1 static char help[] = "Tests 1D discretization tools.\n\n";
2
3 #include <petscdt.h>
4 #include <petscviewer.h>
5 #include <petsc/private/petscimpl.h>
6 #include <petsc/private/dtimpl.h>
7
CheckPoints(const char * name,PetscInt npoints,const PetscReal * points,PetscInt ndegrees,const PetscInt * degrees)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;
12
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 }
32
33 typedef PetscErrorCode (*quadratureFunc)(PetscInt, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal[], PetscReal[]);
34
CheckQuadrature_Basics(PetscInt npoints,PetscReal alpha,PetscReal beta,const PetscReal x[],const PetscReal w[])35 static PetscErrorCode CheckQuadrature_Basics(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal x[], const PetscReal w[])
36 {
37 PetscInt i;
38
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 }
48
CheckQuadrature(PetscInt npoints,PetscReal alpha,PetscReal beta,const PetscReal x[],const PetscReal w[],PetscInt nexact)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;
54
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));
65
66 tol = eps;
67 if (i == j) {
68 PetscReal norm, norm2diff;
69
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;
76
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
81
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);
85
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 }
97
CheckJacobiQuadrature(PetscInt npoints,PetscReal alpha,PetscReal beta,quadratureFunc func,PetscInt nexact)98 static PetscErrorCode CheckJacobiQuadrature(PetscInt npoints, PetscReal alpha, PetscReal beta, quadratureFunc func, PetscInt nexact)
99 {
100 PetscReal *x, *w;
101
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;
114
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;
122
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 }
139
main(int argc,char ** argv)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;
146
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));
156
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));
163
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));
169
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));
182
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;
205
206 for (i = minpoints; i <= maxpoints; i++) {
207 PetscReal a1, b1, a2, b2;
208
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 }
233
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*/
239