Branch data Line data Source code
1 : :
2 : : #include <cmath>
3 : :
4 : : namespace ippl {
5 : : template <typename T, unsigned NumNodes1D, typename ElementType>
6 : 146 : GaussJacobiQuadrature<T, NumNodes1D, ElementType>::GaussJacobiQuadrature(
7 : : const ElementType& ref_element, const T& alpha, const T& beta,
8 : : const size_t& max_newton_iterations, const size_t& min_newton_iterations)
9 : : : Quadrature<T, NumNodes1D, ElementType>(ref_element)
10 : 146 : , alpha_m(alpha)
11 : 146 : , beta_m(beta)
12 : 146 : , max_newton_iterations_m(max_newton_iterations)
13 : 146 : , min_newton_iterations_m(min_newton_iterations) {
14 [ - + ]: 146 : assert(alpha > -1.0 && "alpha >= -1.0 is not satisfied");
15 [ - + ]: 146 : assert(beta > -1.0 && "beta >= -1.0 is not satisfied");
16 [ - + ]: 146 : assert(max_newton_iterations >= 1 && "max_newton_iterations >= 1 is not satisfied");
17 [ - + ]: 146 : assert(min_newton_iterations_m >= 1 && "min_newton_iterations_m >= 1 is not satisfied");
18 [ - + ]: 146 : assert(min_newton_iterations_m <= max_newton_iterations_m
19 : : && "min_newton_iterations_m <= max_newton_iterations_m is not satisfied");
20 : :
21 : 146 : this->degree_m = 2 * NumNodes1D - 1;
22 : :
23 : 146 : this->a_m = -1.0; // start of the domain
24 : 146 : this->b_m = 1.0; // end of the domain
25 : :
26 [ + - ]: 146 : this->integration_nodes_m = Vector<T, NumNodes1D>();
27 [ + - ]: 146 : this->weights_m = Vector<T, NumNodes1D>();
28 : :
29 [ + - ]: 146 : this->computeNodesAndWeights();
30 : 146 : }
31 : :
32 : : template <typename T, unsigned NumNodes1D, typename ElementType>
33 : : typename GaussJacobiQuadrature<T, NumNodes1D, ElementType>::scalar_t
34 : 622 : GaussJacobiQuadrature<T, NumNodes1D, ElementType>::getChebyshevNodes(const size_t& i) const {
35 : 622 : return -Kokkos::cos((2.0 * static_cast<scalar_t>(i) + 1.0) * Kokkos::numbers::pi_v<scalar_t>
36 : 622 : / (2.0 * NumNodes1D));
37 : : }
38 : :
39 : : template <typename T, unsigned NumNodes1D, typename ElementType>
40 : : typename GaussJacobiQuadrature<T, NumNodes1D, ElementType>::scalar_t
41 : 0 : GaussJacobiQuadrature<T, NumNodes1D, ElementType>::getLehrFEMInitialGuess(
42 : : const size_t& i,
43 : : const Vector<GaussJacobiQuadrature<T, NumNodes1D, ElementType>::scalar_t, NumNodes1D>&
44 : : integration_nodes) const {
45 : 0 : const scalar_t alpha = this->alpha_m;
46 : 0 : const scalar_t beta = this->beta_m;
47 : : scalar_t r1;
48 : : scalar_t r2;
49 : : scalar_t r3;
50 [ # # ]: 0 : scalar_t z = (i > 0) ? integration_nodes[i - 1] : 0.0;
51 : :
52 [ # # ]: 0 : if (i == 0) {
53 : : // initial guess for the largest root
54 : 0 : const scalar_t an = alpha / NumNodes1D;
55 : 0 : const scalar_t bn = beta / NumNodes1D;
56 : 0 : r1 = (1.0 + alpha) * (2.78 / (4.0 + NumNodes1D * NumNodes1D) + 0.768 * an / NumNodes1D);
57 : 0 : r2 = 1.0 + 1.48 * an + 0.96 * bn + 0.452 * an * an + 0.83 * an * bn;
58 : 0 : z = 1.0 - r1 / r2;
59 [ # # ]: 0 : } else if (i == 1) {
60 : : // initial guess for the second largest root
61 : 0 : r1 = (4.1 + alpha) / ((1.0 + alpha) * (1.0 + 0.156 * alpha));
62 : 0 : r2 = 1.0 + 0.06 * (NumNodes1D - 8.0) * (1.0 + 0.12 * alpha) / NumNodes1D;
63 : 0 : r3 = 1.0 + 0.012 * beta * (1.0 + 0.25 * Kokkos::abs(alpha)) / NumNodes1D;
64 : 0 : z -= (1.0 - z) * r1 * r2 * r3;
65 [ # # ]: 0 : } else if (i == 2) {
66 : : // initial guess for the third largest root
67 : 0 : r1 = (1.67 + 0.28 * alpha) / (1.0 + 0.37 * alpha);
68 : 0 : r2 = 1.0 + 0.22 * (NumNodes1D - 8.0) / NumNodes1D;
69 : 0 : r3 = 1.0 + 8.0 * beta / ((6.28 + beta) * NumNodes1D * NumNodes1D);
70 : 0 : z -= (integration_nodes[0] - z) * r1 * r2 * r3;
71 [ # # ]: 0 : } else if (i == NumNodes1D - 2) {
72 : : // initial guess for the second smallest root
73 : 0 : r1 = (1.0 + 0.235 * beta) / (0.766 + 0.119 * beta);
74 : 0 : r2 = 1.0 / (1.0 + 0.639 * (NumNodes1D - 4.0) / (1.0 + 0.71 * (NumNodes1D - 4.0)));
75 : 0 : r3 = 1.0 / (1.0 + 20.0 * alpha / ((7.5 + alpha) * NumNodes1D * NumNodes1D));
76 : 0 : z += (z - integration_nodes[NumNodes1D - 4]) * r1 * r2 * r3;
77 [ # # ]: 0 : } else if (i == NumNodes1D - 1) {
78 : : // initial guess for the smallest root
79 : 0 : r1 = (1.0 + 0.37 * beta) / (1.67 + 0.28 * beta);
80 : 0 : r2 = 1.0 / (1.0 + 0.22 * (NumNodes1D - 8.0) / NumNodes1D);
81 : 0 : r3 = 1.0 / (1.0 + 8.0 * alpha / ((6.28 + alpha) * NumNodes1D * NumNodes1D));
82 : 0 : z += (z - integration_nodes[NumNodes1D - 3]) * r1 * r2 * r3;
83 : : } else {
84 : : // initial guess for the other integration_nodes
85 : 0 : z = 3.0 * z - 3.0 * integration_nodes[i - 2] + integration_nodes[i - 3];
86 : : }
87 : 0 : return z;
88 : : }
89 : :
90 : : template <typename T, unsigned NumNodes1D, typename ElementType>
91 : 146 : void GaussJacobiQuadrature<T, NumNodes1D, ElementType>::computeNodesAndWeights() {
92 : : // set the initial guess type
93 : 146 : const InitialGuessType& initial_guess_type = InitialGuessType::Chebyshev;
94 : :
95 : : /**
96 : : * the following algorithm for computing the roots and weights for the Gauss-Jacobi
97 : : * quadrature has been taken from LehrFEM++ (MIT License)
98 : : * https://craffael.github.io/lehrfempp/gauss__quadrature_8cc_source.html
99 : : */
100 : :
101 : 146 : const scalar_t tolerance = 2e-16;
102 : :
103 [ + - ]: 146 : Vector<scalar_t, NumNodes1D> integration_nodes;
104 [ + - ]: 146 : Vector<scalar_t, NumNodes1D> weights;
105 : :
106 : 146 : const scalar_t alpha = this->alpha_m;
107 : 146 : const scalar_t beta = this->beta_m;
108 : :
109 : 146 : const scalar_t alfbet = alpha + beta;
110 : :
111 : : scalar_t a;
112 : : scalar_t b;
113 : : scalar_t c;
114 : : scalar_t p1;
115 : : scalar_t p2;
116 : : scalar_t p3;
117 : : scalar_t pp;
118 : : scalar_t temp;
119 : : scalar_t z;
120 : : scalar_t z1;
121 : :
122 : : // Compute the root of the Jacobi polynomial
123 : :
124 [ + + ]: 712 : for (size_t i = 0; i < NumNodes1D; ++i) {
125 : : // initial guess depending on which root we are computing
126 [ - + ]: 566 : if (initial_guess_type == InitialGuessType::LehrFEM) {
127 [ # # ]: 0 : z = this->getLehrFEMInitialGuess(i, integration_nodes);
128 [ + - ]: 566 : } else if (initial_guess_type == InitialGuessType::Chebyshev) {
129 : 566 : z = -this->getChebyshevNodes(i);
130 : : } else {
131 [ # # # # : 0 : throw IpplException("GaussJacobiQuadrature::computeNodesAndWeights",
# # ]
132 : : "Unknown initial guess type");
133 : : }
134 : :
135 : : // std::cout << NumNodes1D - i - 1 << ", initial guess: " << z << " with "
136 : : // << initial_guess_type << std::endl;
137 : :
138 : 566 : size_t its = 1;
139 : : do {
140 : : // refinement by Newton's method (from LehrFEM++)
141 : 2352 : temp = 2.0 + alfbet;
142 : :
143 : : // Start the recurrence with P_0 and P1 to avoid a division by zero when
144 : : // alpha * beta = 0 or -1
145 : 2352 : p1 = (alpha - beta + temp * z) / 2.0;
146 : 2352 : p2 = 1.0;
147 [ + + ]: 10984 : for (size_t j = 2; j <= NumNodes1D; ++j) {
148 : 8632 : p3 = p2;
149 : 8632 : p2 = p1;
150 : 8632 : temp = 2 * j + alfbet;
151 : 8632 : a = 2 * j * (j + alfbet) * (temp - 2.0);
152 : 8632 : b = (temp - 1.0) * (alpha * alpha - beta * beta + temp * (temp - 2.0) * z);
153 : 8632 : c = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp;
154 : 8632 : p1 = (b * p2 - c * p3) / a;
155 : : }
156 : 2352 : pp = (NumNodes1D * (alpha - beta - temp * z) * p1
157 : 2352 : + 2.0 * (NumNodes1D + alpha) * (NumNodes1D + beta) * p2)
158 : 2352 : / (temp * (1.0 - z * z));
159 : : // p1 is now the desired jacobian polynomial. We next compute pp, its
160 : : // derivative, by a standard relation involving p2, the polynomial of one
161 : : // lower order
162 : 2352 : z1 = z;
163 : 2352 : z = z1 - p1 / pp; // Newtons Formula
164 : :
165 : : // std::cout << "it = " << its << ", i = " << i << ", error: " << Kokkos::abs(z -
166 : : // z1)
167 : : // << std::endl;
168 [ + + + + : 2352 : if (its > this->min_newton_iterations_m && Kokkos::abs(z - z1) <= tolerance) {
+ + ]
169 : 566 : break;
170 : : }
171 : 1786 : ++its;
172 [ + - ]: 1786 : } while (its <= this->max_newton_iterations_m);
173 : :
174 [ - + ]: 566 : if (its > this->max_newton_iterations_m) {
175 : : // TODO switch to inform
176 [ # # # # ]: 0 : std::cout << "Root " << NumNodes1D - i - 1
177 [ # # ]: 0 : << " didn't converge. Tolerance may be too high for data type"
178 [ # # ]: 0 : << std::endl;
179 : : }
180 : :
181 : : // std::cout << "i = " << i << ", result after " << its << " iterations: " << z
182 : : // << std::endl;
183 : :
184 : 566 : integration_nodes[i] = z;
185 : :
186 : : // Compute the weight of the Gauss-Jacobi quadrature
187 : 1132 : weights[i] =
188 : 1132 : Kokkos::exp(Kokkos::lgamma(alpha + NumNodes1D) + Kokkos::lgamma(beta + NumNodes1D)
189 : 566 : - Kokkos::lgamma(NumNodes1D + 1.)
190 : 566 : - Kokkos::lgamma(static_cast<double>(NumNodes1D) + alfbet + 1.0))
191 : 566 : * temp * Kokkos::pow(2.0, alfbet) / (pp * p2);
192 : :
193 : : // store the integration nodes and weights in the correct order
194 : 566 : this->integration_nodes_m[i] = static_cast<T>(-integration_nodes[i]);
195 : 566 : this->weights_m[i] = static_cast<T>(weights[i]);
196 : : }
197 : 146 : }
198 : :
199 : : } // namespace ippl
|
