Line data Source code
1 :
2 : #include <cmath>
3 :
4 : namespace ippl {
5 : template <typename T, unsigned NumNodes1D, typename ElementType>
6 372 : 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 372 : , alpha_m(alpha)
11 372 : , beta_m(beta)
12 372 : , max_newton_iterations_m(max_newton_iterations)
13 372 : , min_newton_iterations_m(min_newton_iterations) {
14 372 : assert(alpha > -1.0 && "alpha >= -1.0 is not satisfied");
15 372 : assert(beta > -1.0 && "beta >= -1.0 is not satisfied");
16 372 : assert(max_newton_iterations >= 1 && "max_newton_iterations >= 1 is not satisfied");
17 372 : assert(min_newton_iterations_m >= 1 && "min_newton_iterations_m >= 1 is not satisfied");
18 372 : assert(min_newton_iterations_m <= max_newton_iterations_m
19 : && "min_newton_iterations_m <= max_newton_iterations_m is not satisfied");
20 :
21 372 : this->degree_m = 2 * NumNodes1D - 1;
22 :
23 372 : this->a_m = -1.0; // start of the domain
24 372 : this->b_m = 1.0; // end of the domain
25 :
26 372 : this->integration_nodes_m = Vector<T, NumNodes1D>();
27 372 : this->weights_m = Vector<T, NumNodes1D>();
28 :
29 372 : this->computeNodesAndWeights();
30 372 : }
31 :
32 : template <typename T, unsigned NumNodes1D, typename ElementType>
33 : typename GaussJacobiQuadrature<T, NumNodes1D, ElementType>::scalar_t
34 1644 : GaussJacobiQuadrature<T, NumNodes1D, ElementType>::getChebyshevNodes(const size_t& i) const {
35 1644 : return -Kokkos::cos((2.0 * static_cast<scalar_t>(i) + 1.0) * Kokkos::numbers::pi_v<scalar_t>
36 1644 : / (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 372 : void GaussJacobiQuadrature<T, NumNodes1D, ElementType>::computeNodesAndWeights() {
92 : // set the initial guess type
93 372 : 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 372 : const scalar_t tolerance = 2e-16;
102 :
103 372 : Vector<scalar_t, NumNodes1D> integration_nodes;
104 372 : Vector<scalar_t, NumNodes1D> weights;
105 :
106 372 : const scalar_t alpha = this->alpha_m;
107 372 : const scalar_t beta = this->beta_m;
108 :
109 372 : 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 1904 : for (size_t i = 0; i < NumNodes1D; ++i) {
125 : // initial guess depending on which root we are computing
126 1532 : if (initial_guess_type == InitialGuessType::LehrFEM) {
127 0 : z = this->getLehrFEMInitialGuess(i, integration_nodes);
128 1532 : } else if (initial_guess_type == InitialGuessType::Chebyshev) {
129 1532 : 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 1532 : size_t its = 1;
139 : do {
140 : // refinement by Newton's method (from LehrFEM++)
141 6624 : 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 6624 : p1 = (alpha - beta + temp * z) / 2.0;
146 6624 : p2 = 1.0;
147 31568 : for (size_t j = 2; j <= NumNodes1D; ++j) {
148 24944 : p3 = p2;
149 24944 : p2 = p1;
150 24944 : temp = 2 * j + alfbet;
151 24944 : a = 2 * j * (j + alfbet) * (temp - 2.0);
152 24944 : b = (temp - 1.0) * (alpha * alpha - beta * beta + temp * (temp - 2.0) * z);
153 24944 : c = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp;
154 24944 : p1 = (b * p2 - c * p3) / a;
155 : }
156 6624 : pp = (NumNodes1D * (alpha - beta - temp * z) * p1
157 6624 : + 2.0 * (NumNodes1D + alpha) * (NumNodes1D + beta) * p2)
158 6624 : / (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 6624 : z1 = z;
163 6624 : z = z1 - p1 / pp; // Newtons Formula
164 :
165 : // std::cout << "it = " << its << ", i = " << i << ", error: " << Kokkos::abs(z -
166 : // z1)
167 : // << std::endl;
168 6624 : if (its > this->min_newton_iterations_m && Kokkos::abs(z - z1) <= tolerance) {
169 1532 : break;
170 : }
171 5092 : ++its;
172 5092 : } while (its <= this->max_newton_iterations_m);
173 :
174 1532 : 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 1532 : integration_nodes[i] = z;
185 :
186 : // Compute the weight of the Gauss-Jacobi quadrature
187 3064 : weights[i] =
188 3064 : Kokkos::exp(Kokkos::lgamma(alpha + NumNodes1D) + Kokkos::lgamma(beta + NumNodes1D)
189 1532 : - Kokkos::lgamma(NumNodes1D + 1.)
190 1532 : - Kokkos::lgamma(static_cast<double>(NumNodes1D) + alfbet + 1.0))
191 1532 : * temp * Kokkos::pow(2.0, alfbet) / (pp * p2);
192 :
193 : // store the integration nodes and weights in the correct order
194 1532 : this->integration_nodes_m[i] = static_cast<T>(-integration_nodes[i]);
195 1532 : this->weights_m[i] = static_cast<T>(weights[i]);
196 : }
197 372 : }
198 :
199 : } // namespace ippl
|
