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, (char *)0, 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*/