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