Actual source code: land_tensors.h
petsc-3.14.6 2021-03-30
1: #define LANDAU_INVSQRT(q) (1./PetscSqrtReal(q))
2: #define LANDAU_SQRT(q) PetscSqrtReal(q)
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: }
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