1 #define LANDAU_INVSQRT(q) (1./PetscSqrtReal(q)) 2 #define LANDAU_SQRT(q) PetscSqrtReal(q) 3 4 /* elliptic functions 5 */ 6 PETSC_DEVICE_FUNC_DECL PetscReal polevl_10(PetscReal x, PetscReal coef[]) 7 { 8 PetscReal ans; 9 PetscInt i; 10 ans = coef[0]; 11 for (i=1; i<11; i++) ans = ans * x + coef[i]; 12 return(ans); 13 } 14 PETSC_DEVICE_FUNC_DECL PetscReal polevl_9(PetscReal x, PetscReal coef[]) 15 { 16 PetscReal ans; 17 PetscInt i; 18 ans = coef[0]; 19 for (i=1; i<10; i++) ans = ans * x + coef[i]; 20 return(ans); 21 } 22 /* 23 * Complete elliptic integral of the second kind 24 */ 25 PETSC_DEVICE_FUNC_DECL void ellipticE(PetscReal x,PetscReal *ret) 26 { 27 #if defined(PETSC_USE_REAL_SINGLE) 28 static PetscReal P2[] = { 29 1.53552577301013293365E-4F, 30 2.50888492163602060990E-3F, 31 8.68786816565889628429E-3F, 32 1.07350949056076193403E-2F, 33 7.77395492516787092951E-3F, 34 7.58395289413514708519E-3F, 35 1.15688436810574127319E-2F, 36 2.18317996015557253103E-2F, 37 5.68051945617860553470E-2F, 38 4.43147180560990850618E-1F, 39 1.00000000000000000299E0F 40 }; 41 static PetscReal Q2[] = { 42 3.27954898576485872656E-5F, 43 1.00962792679356715133E-3F, 44 6.50609489976927491433E-3F, 45 1.68862163993311317300E-2F, 46 2.61769742454493659583E-2F, 47 3.34833904888224918614E-2F, 48 4.27180926518931511717E-2F, 49 5.85936634471101055642E-2F, 50 9.37499997197644278445E-2F, 51 2.49999999999888314361E-1F 52 }; 53 #else 54 static PetscReal P2[] = { 55 1.53552577301013293365E-4, 56 2.50888492163602060990E-3, 57 8.68786816565889628429E-3, 58 1.07350949056076193403E-2, 59 7.77395492516787092951E-3, 60 7.58395289413514708519E-3, 61 1.15688436810574127319E-2, 62 2.18317996015557253103E-2, 63 5.68051945617860553470E-2, 64 4.43147180560990850618E-1, 65 1.00000000000000000299E0 66 }; 67 static PetscReal Q2[] = { 68 3.27954898576485872656E-5, 69 1.00962792679356715133E-3, 70 6.50609489976927491433E-3, 71 1.68862163993311317300E-2, 72 2.61769742454493659583E-2, 73 3.34833904888224918614E-2, 74 4.27180926518931511717E-2, 75 5.85936634471101055642E-2, 76 9.37499997197644278445E-2, 77 2.49999999999888314361E-1 78 }; 79 #endif 80 x = 1 - x; /* where m = 1 - m1 */ 81 *ret = polevl_10(x,P2) - PetscLogReal(x) * (x * polevl_9(x,Q2)); 82 } 83 /* 84 * Complete elliptic integral of the first kind 85 */ 86 PETSC_DEVICE_FUNC_DECL void ellipticK(PetscReal x,PetscReal *ret) 87 { 88 #if defined(PETSC_USE_REAL_SINGLE) 89 static PetscReal P1[] = 90 { 91 1.37982864606273237150E-4F, 92 2.28025724005875567385E-3F, 93 7.97404013220415179367E-3F, 94 9.85821379021226008714E-3F, 95 6.87489687449949877925E-3F, 96 6.18901033637687613229E-3F, 97 8.79078273952743772254E-3F, 98 1.49380448916805252718E-2F, 99 3.08851465246711995998E-2F, 100 9.65735902811690126535E-2F, 101 1.38629436111989062502E0F 102 }; 103 static PetscReal Q1[] = 104 { 105 2.94078955048598507511E-5F, 106 9.14184723865917226571E-4F, 107 5.94058303753167793257E-3F, 108 1.54850516649762399335E-2F, 109 2.39089602715924892727E-2F, 110 3.01204715227604046988E-2F, 111 3.73774314173823228969E-2F, 112 4.88280347570998239232E-2F, 113 7.03124996963957469739E-2F, 114 1.24999999999870820058E-1F, 115 4.99999999999999999821E-1F 116 }; 117 #else 118 static PetscReal P1[] = 119 { 120 1.37982864606273237150E-4, 121 2.28025724005875567385E-3, 122 7.97404013220415179367E-3, 123 9.85821379021226008714E-3, 124 6.87489687449949877925E-3, 125 6.18901033637687613229E-3, 126 8.79078273952743772254E-3, 127 1.49380448916805252718E-2, 128 3.08851465246711995998E-2, 129 9.65735902811690126535E-2, 130 1.38629436111989062502E0 131 }; 132 static PetscReal Q1[] = 133 { 134 2.94078955048598507511E-5, 135 9.14184723865917226571E-4, 136 5.94058303753167793257E-3, 137 1.54850516649762399335E-2, 138 2.39089602715924892727E-2, 139 3.01204715227604046988E-2, 140 3.73774314173823228969E-2, 141 4.88280347570998239232E-2, 142 7.03124996963957469739E-2, 143 1.24999999999870820058E-1, 144 4.99999999999999999821E-1 145 }; 146 #endif 147 x = 1 - x; /* where m = 1 - m1 */ 148 *ret = polevl_10(x,P1) - PetscLogReal(x) * polevl_10(x,Q1); 149 } 150 151 152 /* integration point functions */ 153 /* Evaluates the tensor U=(I-(x-y)(x-y)/(x-y)^2)/|x-y| at point x,y */ 154 /* if x==y we will return zero. This is not the correct result */ 155 /* since the tensor diverges for x==y but when integrated */ 156 /* the divergent part is antisymmetric and vanishes. This is not */ 157 /* trivial, but can be proven. */ 158 #if LANDAU_DIM==3 159 PETSC_DEVICE_FUNC_DECL void LandauTensor3D(const PetscReal x1[], const PetscReal xp, const PetscReal yp, const PetscReal zp, PetscReal U[][3], PetscReal mask) 160 { 161 PetscReal dx[3],inorm3,inorm,inorm2,norm2,x2[] = {xp,yp,zp}; 162 PetscInt d; 163 for (d = 0, norm2 = PETSC_MACHINE_EPSILON; d < 3; ++d) { 164 dx[d] = x2[d] - x1[d]; 165 norm2 += dx[d] * dx[d]; 166 } 167 inorm2 = mask/norm2; 168 inorm = LANDAU_SQRT(inorm2); 169 inorm3 = inorm2*inorm; 170 for (d = 0; d < 3; ++d) U[d][d] = inorm - inorm3 * dx[d] * dx[d]; 171 U[1][0] = U[0][1] = -inorm3 * dx[0] * dx[1]; 172 U[1][2] = U[2][1] = -inorm3 * dx[2] * dx[1]; 173 U[2][0] = U[0][2] = -inorm3 * dx[0] * dx[2]; 174 } 175 #else 176 PETSC_DEVICE_FUNC_DECL void LandauTensor2D(const PetscReal x[], const PetscReal rp, const PetscReal zp, PetscReal Ud[][2], PetscReal Uk[][2], const PetscReal mask) 177 { 178 PetscReal l,s,r=x[0],z=x[1],i1func,i2func,i3func,ks,es,pi4pow,sqrt_1s,r2,rp2,r2prp2,zmzp,zmzp2,tt; 179 //PetscReal mask /* = !!(r!=rp || z!=zp) */; 180 /* !!(zmzp2 > 1.e-12 || (r-rp) > 1.e-12 || (r-rp) < -1.e-12); */ 181 r2=PetscSqr(r); 182 zmzp=z-zp; 183 rp2=PetscSqr(rp); 184 zmzp2=PetscSqr(zmzp); 185 r2prp2=r2+rp2; 186 l = r2 + rp2 + zmzp2; 187 /* if (zmzp2 > PETSC_SMALL) mask = 1; */ 188 /* else if ((tt=(r-rp)) > PETSC_SMALL) mask = 1; */ 189 /* else if (tt < -PETSC_SMALL) mask = 1; */ 190 /* else mask = 0; */ 191 s = mask*2*r*rp/l; /* mask for vectorization */ 192 tt = 1./(1+s); 193 pi4pow = 4*PETSC_PI*LANDAU_INVSQRT(PetscSqr(l)*l); 194 sqrt_1s = LANDAU_SQRT(1.+s); 195 /* sp.ellipe(2.*s/(1.+s)) */ 196 ellipticE(2*s*tt,&es); /* 44 flops * 2 + 75 = 163 flops including 2 logs, 1 sqrt, 1 pow, 21 mult */ 197 /* sp.ellipk(2.*s/(1.+s)) */ 198 ellipticK(2*s*tt,&ks); /* 44 flops + 75 in rest, 21 mult */ 199 /* mask is needed here just for single precision */ 200 i2func = 2./((1-s)*sqrt_1s) * es; 201 i1func = 4./(PetscSqr(s)*sqrt_1s + PETSC_MACHINE_EPSILON) * mask * (ks - (1.+s) * es); 202 i3func = 2./((1-s)*(s)*sqrt_1s + PETSC_MACHINE_EPSILON) * (es - (1-s) * ks); 203 Ud[0][0]= pi4pow*(rp2*i1func+PetscSqr(zmzp)*i2func); 204 Ud[0][1]=Ud[1][0]=Uk[0][1]= -pi4pow*(zmzp)*(r*i2func-rp*i3func); 205 Uk[1][1]=Ud[1][1]= pi4pow*((r2prp2)*i2func-2*r*rp*i3func)*mask; 206 Uk[0][0]= pi4pow*(zmzp2*i3func+r*rp*i1func); 207 Uk[1][0]= -pi4pow*(zmzp)*(r*i3func-rp*i2func); /* 48 mults + 21 + 21 = 90 mults and divs */ 208 } 209 #endif 210