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 8 static PetscErrorCode CheckPoints(const char *name, PetscInt npoints, const PetscReal *points, PetscInt ndegrees, const PetscInt *degrees) { 9 PetscReal *B, *D, *D2; 10 PetscInt i, j; 11 12 PetscFunctionBegin; 13 PetscCall(PetscMalloc3(npoints * ndegrees, &B, npoints * ndegrees, &D, npoints * ndegrees, &D2)); 14 PetscCall(PetscDTLegendreEval(npoints, points, ndegrees, degrees, B, D, D2)); 15 PetscCall(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 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)); 26 } 27 } 28 PetscCall(PetscFree3(B, D, D2)); 29 PetscFunctionReturn(0); 30 } 31 32 typedef PetscErrorCode (*quadratureFunc)(PetscInt, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal[], PetscReal[]); 33 34 static PetscErrorCode CheckQuadrature_Basics(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal x[], const PetscReal w[]) { 35 PetscInt i; 36 37 PetscFunctionBegin; 38 for (i = 1; i < npoints; i++) { 39 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]); 40 } 41 for (i = 0; i < npoints; i++) { 42 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]); 43 } 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode CheckQuadrature(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal x[], const PetscReal w[], PetscInt nexact) { 48 PetscInt i, j, k; 49 PetscReal *Pi, *Pj; 50 PetscReal eps; 51 52 PetscFunctionBegin; 53 eps = PETSC_SMALL; 54 PetscCall(PetscMalloc2(npoints, &Pi, npoints, &Pj)); 55 for (i = 0; i <= nexact; i++) { 56 PetscCall(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 PetscCall(PetscDTJacobiEval(npoints, alpha, beta, x, 1, &j, Pj, NULL, NULL)); 62 63 tol = eps; 64 if (i == j) { 65 PetscReal norm, norm2diff; 66 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; 73 74 PetscCheck((PetscReal)ibeta == beta, PETSC_COMM_SELF, PETSC_ERR_SUP, "lgamma() - math routine is unavailable."); 75 for (k = 0; k < ibeta; k++) I_exact *= (i + 1. + k) / (i + alpha + 1. + k); 76 } 77 #endif 78 79 PetscCall(PetscDTJacobiNorm(alpha, beta, i, &norm)); 80 norm2diff = PetscAbsReal(norm * norm - I_exact); 81 PetscCheck(norm2diff <= eps * I_exact, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Jacobi norm error %g", (double)norm2diff); 82 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 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)); 88 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); 89 } 90 } 91 PetscCall(PetscFree2(Pi, Pj)); 92 PetscFunctionReturn(0); 93 } 94 95 static PetscErrorCode CheckJacobiQuadrature(PetscInt npoints, PetscReal alpha, PetscReal beta, quadratureFunc func, PetscInt nexact) { 96 PetscReal *x, *w; 97 98 PetscFunctionBegin; 99 PetscCall(PetscMalloc2(npoints, &x, npoints, &w)); 100 PetscCall((*func)(npoints, -1., 1., alpha, beta, x, w)); 101 PetscCall(CheckQuadrature_Basics(npoints, alpha, beta, x, w)); 102 PetscCall(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; 110 111 eps = PETSC_SMALL; 112 PetscCall(PetscMalloc2(npoints, &x2, npoints, &w2)); 113 PetscCall((*func)(npoints, -1., 1., alpha, beta, x2, w2)); 114 PetscCall(CheckQuadrature_Basics(npoints, alpha, beta, x2, w2)); 115 PetscCall(CheckQuadrature(npoints, alpha, beta, x2, w2, nexact)); 116 for (i = 0; i < npoints; i++) { 117 PetscReal xdiff, xtol, wdiff, wtol; 118 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 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))); 124 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); 125 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); 126 } 127 PetscCall(PetscFree2(x2, w2)); 128 } 129 /* restore */ 130 PetscDTGaussQuadratureNewton_Internal = (PetscBool)!PetscDTGaussQuadratureNewton_Internal; 131 #endif 132 PetscCall(PetscFree2(x, w)); 133 PetscFunctionReturn(0); 134 } 135 136 int main(int argc, char **argv) { 137 PetscInt degrees[1000], ndegrees, npoints, two; 138 PetscReal points[1000], weights[1000], interval[2]; 139 PetscInt minpoints, maxpoints; 140 PetscBool flg; 141 142 PetscFunctionBeginUser; 143 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 144 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Discretization tools test options", NULL); 145 { 146 ndegrees = 1000; 147 degrees[0] = 0; 148 degrees[1] = 1; 149 degrees[2] = 2; 150 PetscCall(PetscOptionsIntArray("-degrees", "list of degrees to evaluate", "", degrees, &ndegrees, &flg)); 151 152 if (!flg) ndegrees = 3; 153 npoints = 1000; 154 points[0] = 0.0; 155 points[1] = -0.5; 156 points[2] = 1.0; 157 PetscCall(PetscOptionsRealArray("-points", "list of points at which to evaluate", "", points, &npoints, &flg)); 158 159 if (!flg) npoints = 3; 160 two = 2; 161 interval[0] = -1.; 162 interval[1] = 1.; 163 PetscCall(PetscOptionsRealArray("-interval", "interval on which to construct quadrature", "", interval, &two, NULL)); 164 165 minpoints = 1; 166 PetscCall(PetscOptionsInt("-minpoints", "minimum points for thorough Gauss-Jacobi quadrature tests", "", minpoints, &minpoints, NULL)); 167 maxpoints = 30; 168 #if defined(PETSC_USE_REAL_SINGLE) 169 maxpoints = 5; 170 #elif defined(PETSC_USE_REAL___FLOAT128) 171 maxpoints = 20; /* just to make test faster */ 172 #endif 173 PetscCall(PetscOptionsInt("-maxpoints", "maximum points for thorough Gauss-Jacobi quadrature tests", "", maxpoints, &maxpoints, NULL)); 174 } 175 PetscOptionsEnd(); 176 PetscCall(CheckPoints("User-provided points", npoints, points, ndegrees, degrees)); 177 178 PetscCall(PetscDTGaussQuadrature(npoints, interval[0], interval[1], points, weights)); 179 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Quadrature weights\n")); 180 PetscCall(PetscRealView(npoints, weights, PETSC_VIEWER_STDOUT_WORLD)); 181 { 182 PetscReal a = interval[0], b = interval[1], zeroth, first, second; 183 PetscInt i; 184 zeroth = b - a; 185 first = (b * b - a * a) / 2; 186 second = (b * b * b - a * a * a) / 3; 187 for (i = 0; i < npoints; i++) { 188 zeroth -= weights[i]; 189 first -= weights[i] * points[i]; 190 second -= weights[i] * PetscSqr(points[i]); 191 } 192 if (PetscAbs(zeroth) < 1e-10) zeroth = 0.; 193 if (PetscAbs(first) < 1e-10) first = 0.; 194 if (PetscAbs(second) < 1e-10) second = 0.; 195 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Moment error: zeroth=%g, first=%g, second=%g\n", (double)(-zeroth), (double)(-first), (double)(-second))); 196 } 197 PetscCall(CheckPoints("Gauss points", npoints, points, ndegrees, degrees)); 198 { 199 PetscInt i; 200 201 for (i = minpoints; i <= maxpoints; i++) { 202 PetscReal a1, b1, a2, b2; 203 204 #if defined(PETSC_HAVE_LGAMMA) 205 a1 = -0.6; 206 b1 = 1.1; 207 a2 = 2.2; 208 b2 = -0.6; 209 #else 210 a1 = 0.; 211 b1 = 1.; 212 a2 = 2.; 213 b2 = 0.; 214 #endif 215 PetscCall(CheckJacobiQuadrature(i, 0., 0., PetscDTGaussJacobiQuadrature, 2 * i - 1)); 216 PetscCall(CheckJacobiQuadrature(i, a1, b1, PetscDTGaussJacobiQuadrature, 2 * i - 1)); 217 PetscCall(CheckJacobiQuadrature(i, a2, b2, PetscDTGaussJacobiQuadrature, 2 * i - 1)); 218 if (i >= 2) { 219 PetscCall(CheckJacobiQuadrature(i, 0., 0., PetscDTGaussLobattoJacobiQuadrature, 2 * i - 3)); 220 PetscCall(CheckJacobiQuadrature(i, a1, b1, PetscDTGaussLobattoJacobiQuadrature, 2 * i - 3)); 221 PetscCall(CheckJacobiQuadrature(i, a2, b2, PetscDTGaussLobattoJacobiQuadrature, 2 * i - 3)); 222 } 223 } 224 } 225 PetscCall(PetscFinalize()); 226 return 0; 227 } 228 229 /*TEST 230 test: 231 suffix: 1 232 args: -degrees 1,2,3,4,5 -points 0,.2,-.5,.8,.9,1 -interval -.5,1 233 TEST*/ 234