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 { 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 %D 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(0); 31 } 32 33 typedef PetscErrorCode(*quadratureFunc)(PetscInt,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal[],PetscReal[]); 34 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 PetscCheckFalse(x[i] <= x[i-1],PETSC_COMM_SELF,PETSC_ERR_PLIB,"Quadrature points not monotonically increasing, %D points, alpha = %g, beta = %g, i = %D, x[i] = %g, x[i-1] = %g",npoints, (double) alpha, (double) beta, i, x[i], x[i-1]); 42 } 43 for (i = 0; i < npoints; i++) { 44 PetscCheckFalse(w[i] <= 0.,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Quadrature weight not positive, %D points, alpha = %g, beta = %g, i = %D, w[i] = %g",npoints, (double) alpha, (double) beta, i, w[i]); 45 } 46 PetscFunctionReturn(0); 47 } 48 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 PetscCheckFalse((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 PetscCheckFalse(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 %D, alpha %g, beta %g, i %D, j %D, exact %g, err %g\n", npoints, (double) alpha, (double) beta, i, j, (double) I_exact, (double) err)); 91 PetscCheckFalse(err > tol,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrectly integrated P_%D * P_%D using %D 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(0); 96 } 97 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 %D, alpha %g, beta %g, i %D, xdiff/xtol %g, wdiff/wtol %g\n", npoints, (double) alpha, (double) beta, i, (double) xdiff/xtol, (double) wdiff/wtol)); 128 PetscCheckFalse(xdiff > xtol,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatch quadrature point: %D points, alpha = %g, beta = %g, i = %D, xdiff = %g", npoints, (double) alpha, (double) beta, i, (double) xdiff); 129 PetscCheckFalse(wdiff > wtol,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatch quadrature weight: %D points, alpha = %g, beta = %g, i = %D, 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(0); 138 } 139 140 int main(int argc,char **argv) 141 { 142 PetscErrorCode ierr; 143 PetscInt degrees[1000],ndegrees,npoints,two; 144 PetscReal points[1000],weights[1000],interval[2]; 145 PetscInt minpoints, maxpoints; 146 PetscBool flg; 147 148 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 149 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Discretization tools test options",NULL);PetscCall(ierr); 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 ierr = PetscOptionsEnd();PetscCall(ierr); 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