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*/