Actual source code: ex1.c
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: {
10: PetscReal *B, *D, *D2;
11: PetscInt i, j;
13: PetscMalloc3(npoints * ndegrees, &B, npoints * ndegrees, &D, npoints * ndegrees, &D2);
14: PetscDTLegendreEval(npoints, points, ndegrees, degrees, B, D, D2);
15: PetscPrintf(PETSC_COMM_WORLD, "%s\n", name);
16: for (i = 0; i < npoints; i++) {
17: for (j = 0; j < ndegrees; j++) {
18: PetscReal b, d, d2;
19: b = B[i * ndegrees + j];
20: d = D[i * ndegrees + j];
21: d2 = D2[i * ndegrees + j];
22: if (PetscAbsReal(b) < PETSC_SMALL) b = 0;
23: if (PetscAbsReal(d) < PETSC_SMALL) d = 0;
24: if (PetscAbsReal(d2) < PETSC_SMALL) d2 = 0;
25: PetscPrintf(PETSC_COMM_WORLD, "degree %" PetscInt_FMT " at %12.4g: B=%12.4g D=%12.4g D2=%12.4g\n", degrees[j], (double)points[i], (double)b, (double)d, (double)d2);
26: }
27: }
28: PetscFree3(B, D, D2);
29: return 0;
30: }
32: typedef PetscErrorCode (*quadratureFunc)(PetscInt, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal[], PetscReal[]);
34: static PetscErrorCode CheckQuadrature_Basics(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal x[], const PetscReal w[])
35: {
36: PetscInt i;
38: for (i = 1; i < npoints; i++) {
40: }
41: for (i = 0; i < npoints; i++) {
43: }
44: return 0;
45: }
47: static PetscErrorCode CheckQuadrature(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal x[], const PetscReal w[], PetscInt nexact)
48: {
49: PetscInt i, j, k;
50: PetscReal *Pi, *Pj;
51: PetscReal eps;
53: eps = PETSC_SMALL;
54: PetscMalloc2(npoints, &Pi, npoints, &Pj);
55: for (i = 0; i <= nexact; i++) {
56: PetscDTJacobiEval(npoints, alpha, beta, x, 1, &i, Pi, NULL, NULL);
57: for (j = i; j <= nexact - i; j++) {
58: PetscReal I_quad = 0.;
59: PetscReal I_exact = 0.;
60: PetscReal err, tol;
61: PetscDTJacobiEval(npoints, alpha, beta, x, 1, &j, Pj, NULL, NULL);
63: tol = eps;
64: if (i == j) {
65: PetscReal norm, norm2diff;
67: I_exact = PetscPowReal(2.0, alpha + beta + 1.) / (2. * i + alpha + beta + 1.);
68: #if defined(PETSC_HAVE_LGAMMA)
69: I_exact *= PetscExpReal(PetscLGamma(i + alpha + 1.) + PetscLGamma(i + beta + 1.) - (PetscLGamma(i + alpha + beta + 1.) + PetscLGamma(i + 1.)));
70: #else
71: {
72: PetscInt ibeta = (PetscInt)beta;
75: for (k = 0; k < ibeta; k++) I_exact *= (i + 1. + k) / (i + alpha + 1. + k);
76: }
77: #endif
79: PetscDTJacobiNorm(alpha, beta, i, &norm);
80: norm2diff = PetscAbsReal(norm * norm - I_exact);
83: tol = eps * I_exact;
84: }
85: for (k = 0; k < npoints; k++) I_quad += w[k] * (Pi[k] * Pj[k]);
86: err = PetscAbsReal(I_exact - I_quad);
87: PetscInfo(NULL, "npoints %" PetscInt_FMT ", alpha %g, beta %g, i %" PetscInt_FMT ", j %" PetscInt_FMT ", exact %g, err %g\n", npoints, (double)alpha, (double)beta, i, j, (double)I_exact, (double)err);
89: }
90: }
91: PetscFree2(Pi, Pj);
92: return 0;
93: }
95: static PetscErrorCode CheckJacobiQuadrature(PetscInt npoints, PetscReal alpha, PetscReal beta, quadratureFunc func, PetscInt nexact)
96: {
97: 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: PetscInfo(NULL, "npoints %" PetscInt_FMT ", alpha %g, beta %g, i %" PetscInt_FMT ", xdiff/xtol %g, wdiff/wtol %g\n", npoints, (double)alpha, (double)beta, i, (double)(xdiff / xtol), (double)(wdiff / wtol));
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: {
138: PetscInt degrees[1000], ndegrees, npoints, two;
139: PetscReal points[1000], weights[1000], interval[2];
140: PetscInt minpoints, maxpoints;
141: PetscBool flg;
144: PetscInitialize(&argc, &argv, (char *)0, help);
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 0;
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*/