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
|
