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: PetscFunctionBegin;
14: PetscCall(PetscMalloc3(npoints * ndegrees, &B, npoints * ndegrees, &D, npoints * ndegrees, &D2));
15: PetscCall(PetscDTLegendreEval(npoints, points, ndegrees, degrees, B, D, D2));
16: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s\n", name));
17: for (i = 0; i < npoints; i++) {
18: for (j = 0; j < ndegrees; j++) {
19: PetscReal b, d, d2;
20: b = B[i * ndegrees + j];
21: d = D[i * ndegrees + j];
22: d2 = D2[i * ndegrees + j];
23: if (PetscAbsReal(b) < PETSC_SMALL) b = 0;
24: if (PetscAbsReal(d) < PETSC_SMALL) d = 0;
25: if (PetscAbsReal(d2) < PETSC_SMALL) d2 = 0;
26: PetscCall(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));
27: }
28: }
29: PetscCall(PetscFree3(B, D, D2));
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: typedef PetscErrorCode (*quadratureFunc)(PetscInt, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal[], PetscReal[]);
35: static PetscErrorCode CheckQuadrature_Basics(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal x[], const PetscReal w[])
36: {
37: PetscInt i;
39: PetscFunctionBegin;
40: for (i = 1; i < npoints; i++) {
41: PetscCheck(x[i] > x[i - 1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Quadrature points not monotonically increasing, %" PetscInt_FMT " points, alpha = %g, beta = %g, i = %" PetscInt_FMT ", x[i] = %g, x[i-1] = %g", npoints, (double)alpha, (double)beta, i, (double)x[i], (double)x[i - 1]);
42: }
43: for (i = 0; i < npoints; i++) {
44: PetscCheck(w[i] > 0., PETSC_COMM_SELF, PETSC_ERR_PLIB, "Quadrature weight not positive, %" PetscInt_FMT " points, alpha = %g, beta = %g, i = %" PetscInt_FMT ", w[i] = %g", npoints, (double)alpha, (double)beta, i, (double)w[i]);
45: }
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode CheckQuadrature(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal x[], const PetscReal w[], PetscInt nexact)
50: {
51: PetscInt i, j, k;
52: PetscReal *Pi, *Pj;
53: PetscReal eps;
55: PetscFunctionBegin;
56: eps = PETSC_SMALL;
57: PetscCall(PetscMalloc2(npoints, &Pi, npoints, &Pj));
58: for (i = 0; i <= nexact; i++) {
59: PetscCall(PetscDTJacobiEval(npoints, alpha, beta, x, 1, &i, Pi, NULL, NULL));
60: for (j = i; j <= nexact - i; j++) {
61: PetscReal I_quad = 0.;
62: PetscReal I_exact = 0.;
63: PetscReal err, tol;
64: PetscCall(PetscDTJacobiEval(npoints, alpha, beta, x, 1, &j, Pj, NULL, NULL));
66: tol = eps;
67: if (i == j) {
68: PetscReal norm, norm2diff;
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: PetscCheck((PetscReal)ibeta == beta, 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
82: PetscCall(PetscDTJacobiNorm(alpha, beta, i, &norm));
83: norm2diff = PetscAbsReal(norm * norm - I_exact);
84: PetscCheck(norm2diff <= eps * I_exact, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Jacobi norm error %g", (double)norm2diff);
86: tol = eps * I_exact;
87: }
88: for (k = 0; k < npoints; k++) I_quad += w[k] * (Pi[k] * Pj[k]);
89: err = PetscAbsReal(I_exact - I_quad);
90: PetscCall(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));
91: PetscCheck(err <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrectly integrated P_%" PetscInt_FMT " * P_%" PetscInt_FMT " using %" PetscInt_FMT " point rule with alpha = %g, beta = %g: exact %g, err %g", i, j, npoints, (double)alpha, (double)beta, (double)I_exact, (double)err);
92: }
93: }
94: PetscCall(PetscFree2(Pi, Pj));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: static PetscErrorCode CheckJacobiQuadrature(PetscInt npoints, PetscReal alpha, PetscReal beta, quadratureFunc func, PetscInt nexact)
99: {
100: PetscReal *x, *w;
102: PetscFunctionBegin;
103: PetscCall(PetscMalloc2(npoints, &x, npoints, &w));
104: PetscCall((*func)(npoints, -1., 1., alpha, beta, x, w));
105: PetscCall(CheckQuadrature_Basics(npoints, alpha, beta, x, w));
106: PetscCall(CheckQuadrature(npoints, alpha, beta, x, w, nexact));
107: #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
108: /* compare methods of computing quadrature */
109: PetscDTGaussQuadratureNewton_Internal = (PetscBool)!PetscDTGaussQuadratureNewton_Internal;
110: {
111: PetscReal *x2, *w2;
112: PetscReal eps;
113: PetscInt i;
115: eps = PETSC_SMALL;
116: PetscCall(PetscMalloc2(npoints, &x2, npoints, &w2));
117: PetscCall((*func)(npoints, -1., 1., alpha, beta, x2, w2));
118: PetscCall(CheckQuadrature_Basics(npoints, alpha, beta, x2, w2));
119: PetscCall(CheckQuadrature(npoints, alpha, beta, x2, w2, nexact));
120: for (i = 0; i < npoints; i++) {
121: PetscReal xdiff, xtol, wdiff, wtol;
123: xdiff = PetscAbsReal(x[i] - x2[i]);
124: wdiff = PetscAbsReal(w[i] - w2[i]);
125: xtol = eps * (1. + PetscMin(PetscAbsReal(x[i]), 1. - PetscAbsReal(x[i])));
126: wtol = eps * (1. + w[i]);
127: PetscCall(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)));
128: PetscCheck(xdiff <= xtol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch quadrature point: %" PetscInt_FMT " points, alpha = %g, beta = %g, i = %" PetscInt_FMT ", xdiff = %g", npoints, (double)alpha, (double)beta, i, (double)xdiff);
129: PetscCheck(wdiff <= wtol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatch quadrature weight: %" PetscInt_FMT " points, alpha = %g, beta = %g, i = %" PetscInt_FMT ", wdiff = %g", npoints, (double)alpha, (double)beta, i, (double)wdiff);
130: }
131: PetscCall(PetscFree2(x2, w2));
132: }
133: /* restore */
134: PetscDTGaussQuadratureNewton_Internal = (PetscBool)!PetscDTGaussQuadratureNewton_Internal;
135: #endif
136: PetscCall(PetscFree2(x, w));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: int main(int argc, char **argv)
141: {
142: PetscInt degrees[1000], ndegrees, npoints, two;
143: PetscReal points[1000], weights[1000], interval[2];
144: PetscInt minpoints, maxpoints;
145: PetscBool flg;
147: PetscFunctionBeginUser;
148: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
149: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Discretization tools test options", NULL);
150: {
151: ndegrees = 1000;
152: degrees[0] = 0;
153: degrees[1] = 1;
154: degrees[2] = 2;
155: PetscCall(PetscOptionsIntArray("-degrees", "list of degrees to evaluate", "", degrees, &ndegrees, &flg));
157: if (!flg) ndegrees = 3;
158: npoints = 1000;
159: points[0] = 0.0;
160: points[1] = -0.5;
161: points[2] = 1.0;
162: PetscCall(PetscOptionsRealArray("-points", "list of points at which to evaluate", "", points, &npoints, &flg));
164: if (!flg) npoints = 3;
165: two = 2;
166: interval[0] = -1.;
167: interval[1] = 1.;
168: PetscCall(PetscOptionsRealArray("-interval", "interval on which to construct quadrature", "", interval, &two, NULL));
170: minpoints = 1;
171: PetscCall(PetscOptionsInt("-minpoints", "minimum points for thorough Gauss-Jacobi quadrature tests", "", minpoints, &minpoints, NULL));
172: maxpoints = 30;
173: #if defined(PETSC_USE_REAL_SINGLE)
174: maxpoints = 5;
175: #elif defined(PETSC_USE_REAL___FLOAT128)
176: maxpoints = 20; /* just to make test faster */
177: #endif
178: PetscCall(PetscOptionsInt("-maxpoints", "maximum points for thorough Gauss-Jacobi quadrature tests", "", maxpoints, &maxpoints, NULL));
179: }
180: PetscOptionsEnd();
181: PetscCall(CheckPoints("User-provided points", npoints, points, ndegrees, degrees));
183: PetscCall(PetscDTGaussQuadrature(npoints, interval[0], interval[1], points, weights));
184: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Quadrature weights\n"));
185: PetscCall(PetscRealView(npoints, weights, PETSC_VIEWER_STDOUT_WORLD));
186: {
187: PetscReal a = interval[0], b = interval[1], zeroth, first, second;
188: PetscInt i;
189: zeroth = b - a;
190: first = (b * b - a * a) / 2;
191: second = (b * b * b - a * a * a) / 3;
192: for (i = 0; i < npoints; i++) {
193: zeroth -= weights[i];
194: first -= weights[i] * points[i];
195: second -= weights[i] * PetscSqr(points[i]);
196: }
197: if (PetscAbs(zeroth) < 1e-10) zeroth = 0.;
198: if (PetscAbs(first) < 1e-10) first = 0.;
199: if (PetscAbs(second) < 1e-10) second = 0.;
200: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Moment error: zeroth=%g, first=%g, second=%g\n", (double)(-zeroth), (double)(-first), (double)(-second)));
201: }
202: PetscCall(CheckPoints("Gauss points", npoints, points, ndegrees, degrees));
203: {
204: PetscInt i;
206: for (i = minpoints; i <= maxpoints; i++) {
207: PetscReal a1, b1, a2, b2;
209: #if defined(PETSC_HAVE_LGAMMA)
210: a1 = -0.6;
211: b1 = 1.1;
212: a2 = 2.2;
213: b2 = -0.6;
214: #else
215: a1 = 0.;
216: b1 = 1.;
217: a2 = 2.;
218: b2 = 0.;
219: #endif
220: PetscCall(CheckJacobiQuadrature(i, 0., 0., PetscDTGaussJacobiQuadrature, 2 * i - 1));
221: PetscCall(CheckJacobiQuadrature(i, a1, b1, PetscDTGaussJacobiQuadrature, 2 * i - 1));
222: PetscCall(CheckJacobiQuadrature(i, a2, b2, PetscDTGaussJacobiQuadrature, 2 * i - 1));
223: if (i >= 2) {
224: PetscCall(CheckJacobiQuadrature(i, 0., 0., PetscDTGaussLobattoJacobiQuadrature, 2 * i - 3));
225: PetscCall(CheckJacobiQuadrature(i, a1, b1, PetscDTGaussLobattoJacobiQuadrature, 2 * i - 3));
226: PetscCall(CheckJacobiQuadrature(i, a2, b2, PetscDTGaussLobattoJacobiQuadrature, 2 * i - 3));
227: }
228: }
229: }
230: PetscCall(PetscFinalize());
231: return 0;
232: }
234: /*TEST
235: test:
236: suffix: 1
237: args: -degrees 1,2,3,4,5 -points 0,.2,-.5,.8,.9,1 -interval -.5,1
238: TEST*/