DiFfRG
Loading...
Searching...
No Matches
loop_integrals.hh
Go to the documentation of this file.
1#pragma once
2
3// standard library
4#include <cmath>
5
6// external libraries
7#include <deal.II/base/quadrature_lib.h>
8#include <tbb/tbb.h>
9
10// DiFfRG
11#include <DiFfRG/common/utils.hh>
12
13namespace DiFfRG
14{
15 using namespace dealii;
16
17 namespace LoopIntegrals
18 {
34 template <typename NT, int d, typename FUN>
35 NT integrate(const FUN &fun, const QGauss<1> &x_quadrature, const double x_extent, const double k)
36 {
37 const auto &x_q_p = x_quadrature.get_points();
38 const auto &x_q_w = x_quadrature.get_weights();
39 const int x_size = x_quadrature.size();
40
41 static_assert(d > 0, "Dimension must be greater than zero");
42 // Dimension of the spatial integral
43 const double S_d = 2. * std::pow(M_PI, d / 2.) / std::tgammal(d / 2.);
44 // Prefactor for the spatial integral
45 const double prefactor = S_d // angular integral
46 * powr<-d>(2. * M_PI); // fourier factors
47
48 const NT result = tbb::parallel_reduce(
49 tbb::blocked_range<int>(0, x_size), NT(0),
50 [&](const tbb::blocked_range<int> &r, NT running_total) -> NT {
51 for (int x_it = r.begin(); x_it < r.end(); x_it++) {
52 const double x = x_q_p[x_it][0] * x_extent;
53 const double x_weight = x_q_w[x_it] * x_extent;
54 const double q2 = x * powr<2>(k);
55
56 running_total +=
57 x_weight * prefactor *
58 (0.5 * powr<d>(k) * std::pow(x, (d - 2) / 2.)) // integral over x = p^2 / k^2 in d dimensions
59 * fun(q2); // integrand
60 }
61 return running_total;
62 },
63 [](const NT a, const NT b) { return a + b; });
64 return result;
65 }
66
83 template <typename NT, int d, typename FUN>
84 NT angle_integrate(const FUN &fun, const QGauss<1> &x_quadrature, const double x_extent, const double k,
85 const QGauss<1> &cos_quadrature)
86 {
87 static_assert(d > 1, "Dimension must be greater than one, otherwise the angular integral not defined.");
88
89 const auto &x_q_p = x_quadrature.get_points();
90 const auto &x_q_w = x_quadrature.get_weights();
91 const int x_size = x_quadrature.size();
92 const auto &cos_q_p = cos_quadrature.get_points();
93 const auto &cos_q_w = cos_quadrature.get_weights();
94 const int cos_size = cos_quadrature.size();
95
96 static_assert(d > 1, "Dimension must be greater than one");
97 // Area of unit sphere in d dimensions
98 const double S_d = 2. * std::pow(M_PI, d / 2.) / std::tgamma(d / 2.);
99 // Prefactor for the spatial integral
100 const double prefactor = S_d // angular integral
101 * powr<-d>(2. * M_PI); // fourier factors
102
103 const NT result = tbb::parallel_reduce(
104 tbb::blocked_range2d<int>(0, x_size, 0, cos_size), NT(0),
105 [&](const tbb::blocked_range2d<int> &r, NT running_total) -> NT {
106 for (int x_it = r.rows().begin(); x_it < r.rows().end(); x_it++) {
107 const double x = x_q_p[x_it][0] * x_extent;
108 const double x_weight = x_q_w[x_it] * x_extent;
109 const double q2 = x * powr<2>(k);
110
111 for (int cos_it = r.cols().begin(); cos_it < r.cols().end(); cos_it++) {
112 const double cos = (cos_q_p[cos_it][0] - 0.5) * 2.;
113 const double cos_weight = cos_q_w[cos_it] * 2.;
114
115 running_total +=
116 x_weight * prefactor *
117 (0.5 * powr<d>(k) * std::pow(x, (d - 2) / 2.)) // integral over x = p^2 / k^2 in d dimensions
118 * 0.5 * cos_weight // integral over cos(theta), 0.5 removes the factor from the angular integral
119 * fun(q2, cos); // integrand
120 }
121 }
122 return running_total;
123 },
124 [](const NT a, const NT b) { return a + b; });
125 return result;
126 }
127
144 template <typename NT, int d, typename FUN>
145 NT two_angle_integrate(const FUN &fun, const QGauss<1> &x_quadrature, const double x_extent, const double k,
146 const QGauss<1> &cos_quadrature)
147 {
148 static_assert(d > 2, "Dimension must be greater than 2, otherwise the angular integral not defined.");
149
150 const auto &x_q_p = x_quadrature.get_points();
151 const auto &x_q_w = x_quadrature.get_weights();
152 const int x_size = x_quadrature.size();
153 const auto &cos_q_p = cos_quadrature.get_points();
154 const auto &cos_q_w = cos_quadrature.get_weights();
155 const int cos_size = cos_quadrature.size();
156 const auto &phi_q_p = cos_quadrature.get_points();
157 const auto &phi_q_w = cos_quadrature.get_weights();
158 const int phi_size = cos_quadrature.size();
159
160 static_assert(d > 1, "Dimension must be greater than one");
161 // Area of unit sphere in d dimensions
162 const double S_d = 2. * std::pow(M_PI, d / 2.) / std::tgamma(d / 2.);
163 // Prefactor for the spatial integral
164 const double prefactor = S_d // angular integral
165 * powr<-d>(2. * M_PI); // fourier factors
166
167 const NT result = tbb::parallel_reduce(
168 tbb::blocked_range3d<int>(0, x_size, 0, cos_size, 0, phi_size), NT(0),
169 [&](const tbb::blocked_range3d<int> &r, NT running_total) -> NT {
170 for (int x_it = r.pages().begin(); x_it < r.pages().end(); x_it++) {
171 const double x = x_q_p[x_it][0] * x_extent;
172 const double x_weight = x_q_w[x_it] * x_extent;
173 const double q2 = x * powr<2>(k);
174
175 for (int cos_it = r.rows().begin(); cos_it < r.rows().end(); cos_it++) {
176 const double cos = (cos_q_p[cos_it][0] - 0.5) * 2.;
177 const double cos_weight = cos_q_w[cos_it] * 2.;
178
179 for (int phi_it = r.cols().begin(); phi_it < r.cols().end(); phi_it++) {
180 const double phi = phi_q_p[phi_it][0] * 2. * M_PI;
181 const double phi_weight = cos_q_w[cos_it] * 2. * M_PI;
182
183 running_total +=
184 x_weight * prefactor *
185 (0.5 * powr<d>(k) * std::pow(x, (d - 2) / 2.)) // integral over x = p^2 / k^2 in d dimensions
186 * 0.5 * cos_weight // integral over cos(theta), 0.5 removes the factor from the angular integral
187 * phi_weight / (2. * M_PI) // integral over phi, 2pi removes the factor from the angular integral
188 * fun(q2, cos, phi); // integrand
189 }
190 }
191 }
192 return running_total;
193 },
194 [](const NT a, const NT b) { return a + b; });
195 return result;
196 }
197
214 template <typename NT, int d, typename FUN>
215 NT three_angle_integrate(const FUN &fun, const QGauss<1> &x_quadrature, const double x_extent, const double k,
216 const QGauss<1> &cos_quadrature)
217 {
218 static_assert(d > 3, "Dimension must be greater than 3, otherwise the angular integral not defined.");
219
220 const auto &x_q_p = x_quadrature.get_points();
221 const auto &x_q_w = x_quadrature.get_weights();
222 const int x_size = x_quadrature.size();
223 const auto &cos_q_p = cos_quadrature.get_points();
224 const auto &cos_q_w = cos_quadrature.get_weights();
225 const int cos_size = cos_quadrature.size();
226 const auto &phi_q_p = cos_quadrature.get_points();
227 const auto &phi_q_w = cos_quadrature.get_weights();
228 const int phi_size = cos_quadrature.size();
229
230 static_assert(d > 1, "Dimension must be greater than one");
231 // Area of unit sphere in d dimensions
232 const double S_d = 2. * std::pow(M_PI, d / 2.) / std::tgamma(d / 2.);
233 const double S_4 = 2. * std::pow(M_PI, 4. / 2.) / std::tgamma(4. / 2.);
234 // Prefactor for the spatial integral
235 const double prefactor = S_d // angular integral
236 * powr<-d>(2. * M_PI); // fourier factors
237
238 const NT result = tbb::parallel_reduce(
239 tbb::blocked_range3d<int>(0, x_size, 0, cos_size, 0, cos_size), NT(0),
240 [&](const tbb::blocked_range3d<int> &r, NT running_total) -> NT {
241 for (int x_it = r.pages().begin(); x_it < r.pages().end(); x_it++) {
242 const double x = x_q_p[x_it][0] * x_extent;
243 const double x_weight = x_q_w[x_it] * x_extent;
244 const double q2 = x * powr<2>(k);
245
246 for (int cos1_it = r.rows().begin(); cos1_it < r.rows().end(); cos1_it++) {
247 const double cos1 = (cos_q_p[cos1_it][0] - 0.5) * 2.;
248 const double cos1_weight =
249 cos_q_w[cos1_it] * 2. * std::sqrt(1. - powr<2>(cos1)); // integration jacobian
250
251 for (int cos2_it = r.cols().begin(); cos2_it < r.cols().end(); cos2_it++) {
252 const double cos2 = (cos_q_p[cos2_it][0] - 0.5) * 2.;
253 const double cos2_weight = cos_q_w[cos2_it] * 2.;
254
255 for (int phi_it = 0; phi_it < phi_size; phi_it++) {
256 const double phi = phi_q_p[phi_it][0] * 2. * M_PI;
257 const double phi_weight = phi_q_w[phi_it] * 2. * M_PI;
258
259 running_total +=
260 x_weight * prefactor *
261 (0.5 * powr<d>(k) * std::pow(x, (d - 2) / 2.)) // integral over x = p^2 / k^2 in d dimensions
262 * cos1_weight // integral over cos(theta1)
263 * cos2_weight // integral over cos(theta2)
264 * phi_weight /
265 S_4 // integral over phi, S_3 = 2 pi^2 removes the factor from the angular integral
266 * fun(q2, cos1, cos2, phi); // integrand
267 }
268 }
269 }
270 }
271 return running_total;
272 },
273 [](const NT a, const NT b) { return a + b; });
274 return result;
275 }
276
295 template <typename NT, int d, typename FUN>
296 NT spatial_integrate_and_integrate(const FUN &fun, const QGauss<1> &x_quadrature, const double x_extent,
297 const double k, const QGauss<1> &m_quadrature, const double m_extent)
298 {
299 const auto &x_q_p = x_quadrature.get_points();
300 const auto &x_q_w = x_quadrature.get_weights();
301 const int x_size = x_quadrature.size();
302 const auto &m_q_p = m_quadrature.get_points();
303 const auto &m_q_w = m_quadrature.get_weights();
304 const int m_size = m_quadrature.size();
305
306 static_assert(d > 0, "Dimension must be greater than zero");
307 // Dimension of the spatial integral
308 constexpr int ds = d - 1;
309 // Area of unit sphere in ds=d-1 dimensions
310 const double S_ds = 2. * std::pow(M_PI, ds / 2.) / std::tgamma(ds / 2.);
311 // Prefactor for the spatial integral
312 const double prefactor = S_ds // angular integral
313 * powr<-ds>(2. * M_PI); // fourier factors
314
315 const NT result = tbb::parallel_reduce(
316 tbb::blocked_range2d<int>(0, x_size, 0, m_size), NT(0),
317 [&](const tbb::blocked_range2d<int> &r, NT running_total) -> NT {
318 for (int x_it = r.rows().begin(); x_it < r.rows().end(); x_it++) {
319 const double x = x_q_p[x_it][0] * x_extent;
320 const double x_weight = x_q_w[x_it] * x_extent;
321 const double q2 = x * powr<2>(k);
322
323 // Do a "matsubara" (rather p0)-integral
324 for (int q0_it = r.cols().begin(); q0_it < r.cols().end(); q0_it++) {
325 const double q0 = m_q_p[q0_it][0] * m_extent;
326 const double m_weight = m_q_w[q0_it] * m_extent;
327
328 running_total += (m_weight / (2. * M_PI)) // Matsubara integral weight
329 * x_weight * prefactor *
330 (0.5 * powr<ds>(k) * std::pow(x, (ds - 2) / 2.)) // integral over x = p^2 / k^2
331 * (fun(q2, q0) + fun(q2, -q0)); // integrand
332 }
333 }
334 return running_total;
335 },
336 [](const NT a, const NT b) { return a + b; });
337 return result;
338 }
339
359 template <typename NT, int d, typename FUN>
360 NT spatial_angle_integrate_and_integrate(const FUN &fun, const QGauss<1> &x_quadrature, const double x_extent,
361 const double k, const QGauss<1> &cos_quadrature,
362 const QGauss<1> &m_quadrature, const double m_extent)
363 {
364 const auto &x_q_p = x_quadrature.get_points();
365 const auto &x_q_w = x_quadrature.get_weights();
366 const int x_size = x_quadrature.size();
367 const auto &cos_q_p = cos_quadrature.get_points();
368 const auto &cos_q_w = cos_quadrature.get_weights();
369 const int cos_size = cos_quadrature.size();
370 const auto &m_q_p = m_quadrature.get_points();
371 const auto &m_q_w = m_quadrature.get_weights();
372 const int m_size = m_quadrature.size();
373
374 static_assert(d > 0, "Dimension must be greater than zero");
375 // Dimension of the spatial integral
376 constexpr int ds = d - 1;
377 // Area of unit sphere in ds=d-1 dimensions
378 const double S_ds = 2. * std::pow(M_PI, ds / 2.) / std::tgamma(ds / 2.);
379 // Prefactor for the spatial integral
380 const double prefactor = S_ds // angular integral
381 * powr<-ds>(2. * M_PI); // fourier factors
382
383 const NT result = tbb::parallel_reduce(
384 tbb::blocked_range3d<int>(0, x_size, 0, cos_size, 0, m_size), NT(0),
385 [&](const tbb::blocked_range3d<int> &r, NT running_total) -> NT {
386 for (int x_it = r.pages().begin(); x_it < r.pages().end(); x_it++) {
387 const double x = x_q_p[x_it][0] * x_extent;
388 const double x_weight = x_q_w[x_it] * x_extent;
389 const double q2 = x * powr<2>(k);
390
391 for (int cos_it = r.rows().begin(); cos_it < r.rows().end(); cos_it++) {
392 const double cos = (cos_q_p[cos_it][0] - 0.5) * 2.;
393 const double cos_weight = cos_q_w[cos_it] * 2.;
394
395 // Do a "matsubara" (rather q0)-integral
396 for (int q0_it = r.cols().begin(); q0_it < r.cols().end(); q0_it++) {
397 const double q0 = m_q_p[q0_it][0] * m_extent;
398 const double m_weight = m_q_w[q0_it] * m_extent;
399
400 running_total +=
401 (m_weight / (2. * M_PI)) // Matsubara integral weight
402 * x_weight * prefactor *
403 (0.5 * powr<ds>(k) * std::pow(x, (ds - 2) / 2.)) // integral over x = p^2 / k^2
404 * 0.5 * cos_weight // integral over cos(theta), 0.5 removes the factor from the angular integral
405 * (fun(q2, cos, q0) + fun(q2, cos, -q0)); // integrand
406 }
407 }
408 }
409 return running_total;
410 },
411 [](const NT a, const NT b) { return a + b; });
412 return result;
413 }
414
435 template <typename NT, int d, typename FUN>
436 NT spatial_integrate_and_sum(const FUN &fun, const QGauss<1> &x_quadrature, const double x_extent, const double k,
437 const int m_order, const QGauss<1> &m_quadrature, const double m_extent,
438 const double T)
439 {
440 if (is_close(T, 0.))
441 return spatial_integrate_and_integrate<NT, d, FUN>(fun, x_quadrature, x_extent, k, m_quadrature, m_extent);
442
443 const auto &x_q_p = x_quadrature.get_points();
444 const auto &x_q_w = x_quadrature.get_weights();
445 const int x_size = x_quadrature.size();
446 const auto &m_q_p = m_quadrature.get_points();
447 const auto &m_q_w = m_quadrature.get_weights();
448 const int m_size = m_quadrature.size();
449
450 static_assert(d > 0, "Dimension must be greater than zero");
451 // Dimension of the spatial integral
452 constexpr int ds = d - 1;
453 // Area of unit sphere in ds=d-1 dimensions
454 const double S_ds = 2. * std::pow(M_PI, ds / 2.) / std::tgamma(ds / 2.);
455 // Prefactor for the spatial integral
456 const double prefactor = S_ds // angular integral
457 * powr<-ds>(2. * M_PI); // fourier factors
458
459 // Do a matsubara sum for the first m_order summands
460 const NT result_sum = tbb::parallel_reduce(
461 tbb::blocked_range2d<int>(0, x_size, -m_order, m_order + 1), NT(0),
462 [&](const tbb::blocked_range2d<int> &r, NT running_total) -> NT {
463 for (int x_it = r.rows().begin(); x_it < r.rows().end(); x_it++) {
464 const double x = x_q_p[x_it][0] * x_extent;
465 const double x_weight = x_q_w[x_it] * x_extent;
466 const double q2 = x * powr<2>(k);
467
468 for (int q0_it = r.cols().begin(); q0_it < r.cols().end(); q0_it++) {
469 const double q0 = 2. * M_PI * T * q0_it;
470
471 running_total += T // Matsubara sum weight
472 * x_weight * prefactor *
473 (0.5 * powr<ds>(k) * std::pow(x, (ds - 2) / 2.)) // integral over x = p^2 / k^2
474 * fun(q2, q0); // integrand
475 }
476 }
477 return running_total;
478 },
479 [](const NT a, const NT b) { return a + b; });
480
481 // Approximate the rest of the sum by an integral
482 const NT result_integral = tbb::parallel_reduce(
483 tbb::blocked_range2d<int>(0, x_size, 0, m_size), NT(0),
484 [&](const tbb::blocked_range2d<int> &r, NT running_total) -> NT {
485 for (int x_it = r.rows().begin(); x_it < r.rows().end(); x_it++) {
486 const double x = x_q_p[x_it][0] * x_extent;
487 const double x_weight = x_q_w[x_it] * x_extent;
488 const double q2 = x * powr<2>(k);
489
490 for (int q0_it = r.cols().begin(); q0_it < r.cols().end(); q0_it++) {
491 const double q0 =
492 (2. * m_order * M_PI * T) + m_q_p[q0_it][0] * (m_extent - 2. * (m_order + 1) * M_PI * T);
493 const double m_weight = m_q_w[q0_it] * (m_extent - 2. * m_order * M_PI * T);
494
495 running_total += (1. * m_weight / (2. * M_PI)) // Matsubara integral weight
496 * x_weight * prefactor *
497 (0.5 * powr<ds>(k) * std::pow(x, (ds - 2) / 2.)) // integral over x = p^2 / k^2
498 * (fun(q2, q0) + fun(q2, -q0)); // integrand
499 }
500 }
501 return running_total;
502 },
503 [](const NT a, const NT b) { return a + b; });
504 return result_sum + result_integral;
505 } // namespace LoopIntegrals
506
528 template <typename NT, int d, typename FUN>
529 NT spatial_angle_integrate_and_sum(const FUN &fun, const QGauss<1> &x_quadrature, const double x_extent,
530 const double k, const QGauss<1> &cos_quadrature, const int m_order,
531 const QGauss<1> &m_quadrature, const double m_extent, const double T)
532 {
533 if (is_close(T, 0.))
534 return spatial_angle_integrate_and_integrate<NT, d, FUN>(fun, x_quadrature, x_extent, k, cos_quadrature,
535 m_quadrature, m_extent);
536
537 const auto &x_q_p = x_quadrature.get_points();
538 const auto &x_q_w = x_quadrature.get_weights();
539 const int x_size = x_quadrature.size();
540 const auto &cos_q_p = cos_quadrature.get_points();
541 const auto &cos_q_w = cos_quadrature.get_weights();
542 const int cos_size = cos_quadrature.size();
543 const auto &m_q_p = m_quadrature.get_points();
544 const auto &m_q_w = m_quadrature.get_weights();
545 const int m_size = m_quadrature.size();
546
547 static_assert(d > 0, "Dimension must be greater than zero");
548 // Dimension of the spatial integral
549 constexpr int ds = d - 1;
550 // Area of unit sphere in ds=d-1 dimensions
551 const double S_ds = 2. * std::pow(M_PI, ds / 2.) / std::tgamma(ds / 2.);
552 // Prefactor for the spatial integral
553 const double prefactor = S_ds // angular integral
554 * powr<-ds>(2. * M_PI); // fourier factors
555
556 // Do a matsubara sum for the first m_order summands
557 const NT result_sum = tbb::parallel_reduce(
558 tbb::blocked_range3d<int>(0, x_size, 0, cos_size, -m_order, m_order + 1), NT(0),
559 [&](const tbb::blocked_range3d<int> &r, NT running_total) -> NT {
560 for (int x_it = r.pages().begin(); x_it < r.pages().end(); x_it++) {
561 const double x = x_q_p[x_it][0] * x_extent;
562 const double x_weight = x_q_w[x_it] * x_extent;
563 const double q2 = x * powr<2>(k);
564
565 for (int cos_it = r.rows().begin(); cos_it < r.rows().end(); cos_it++) {
566 const double cos = (cos_q_p[cos_it][0] - 0.5) * 2.;
567 const double cos_weight = cos_q_w[cos_it] * 2.;
568
569 for (int q0_it = r.cols().begin(); q0_it < r.cols().end(); q0_it++) {
570 const double q0 = 2. * M_PI * T * q0_it;
571
572 running_total +=
573 T // Matsubara sum weight
574 * x_weight * prefactor *
575 (0.5 * powr<ds>(k) * std::pow(x, (ds - 2) / 2.)) // integral over x = p^2 / k^2
576 * 0.5 * cos_weight // integral over cos(theta), 0.5 removes the factor from the angular integral
577 * fun(q2, cos, q0); // integrand
578 }
579 }
580 }
581 return running_total;
582 },
583 [](const NT a, const NT b) { return a + b; });
584
585 // Approximate the rest of the sum by an integral
586 const NT result_integral = tbb::parallel_reduce(
587 tbb::blocked_range3d<int>(0, x_size, 0, cos_size, 0, m_size), NT(0),
588 [&](const tbb::blocked_range3d<int> &r, NT running_total) -> NT {
589 for (int x_it = r.pages().begin(); x_it < r.pages().end(); x_it++) {
590 const double x = x_q_p[x_it][0] * x_extent;
591 const double x_weight = x_q_w[x_it] * x_extent;
592 const double q2 = x * powr<2>(k);
593
594 for (int cos_it = r.rows().begin(); cos_it < r.rows().end(); cos_it++) {
595 const double cos = (cos_q_p[cos_it][0] - 0.5) * 2.;
596 const double cos_weight = cos_q_w[cos_it] * 2.;
597
598 for (int q0_it = r.cols().begin(); q0_it < r.cols().end(); q0_it++) {
599 const double q0 =
600 (2. * m_order * M_PI * T) + m_q_p[q0_it][0] * (m_extent - 2. * (m_order + 1) * M_PI * T);
601 const double m_weight = m_q_w[q0_it] * (m_extent - 2. * m_order * M_PI * T);
602
603 running_total +=
604 (1. * m_weight / (2. * M_PI)) // Matsubara integral weight
605 * x_weight * prefactor *
606 (0.5 * powr<ds>(k) * std::pow(x, (ds - 2) / 2.)) // integral over x = p^2 / k^2
607 * 0.5 * cos_weight // integral over cos(theta), 0.5 removes the factor from the angular integral
608 * (fun(q2, cos, q0) + fun(q2, cos, -q0)); // integrand
609 }
610 }
611 }
612 return running_total;
613 },
614 [](const NT a, const NT b) { return a + b; });
615 return result_sum + result_integral;
616 }
617
618 template <typename NT, typename FUN>
619 NT sum(const FUN &fun, const int m_order, const QGauss<1> &m_quadrature, const double m_extent, const double T)
620 {
621 const auto &m_q_p = m_quadrature.get_points();
622 const auto &m_q_w = m_quadrature.get_weights();
623 const int m_size = m_quadrature.size();
624
625 // Do a matsubara sum for the first m_order summands
626 const NT result_sum = is_close(T, 0.) ? 0.
627 : tbb::parallel_reduce(
628 tbb::blocked_range<int>(-m_order, m_order + 1), NT(0),
629 [&](const tbb::blocked_range<int> &r, NT running_total) -> NT {
630 for (int q0_it = r.begin(); q0_it < r.end(); q0_it++) {
631 const double q0 = 2. * M_PI * T * q0_it;
632 running_total += T * fun(q0);
633 }
634 return running_total;
635 },
636 [](const NT a, const NT b) { return a + b; });
637
638 // Approximate the rest of the sum by an integral
639 const NT result_integral = tbb::parallel_reduce(
640 tbb::blocked_range<int>(0, m_size), NT(0),
641 [&](const tbb::blocked_range<int> &r, NT running_total) -> NT {
642 for (int q0_it = r.begin(); q0_it < r.end(); q0_it++) {
643 const double q0 =
644 (2. * m_order * M_PI * T) + m_q_p[q0_it][0] * (m_extent - 2. * (m_order + 1) * M_PI * T);
645 const double m_weight = m_q_w[q0_it] * (m_extent - 2. * m_order * M_PI * T);
646
647 running_total += (1. * m_weight / (2. * M_PI)) // Matsubara integral weight
648 * (fun(q0) + fun(-q0)); // integrand
649 }
650 return running_total;
651 },
652 [](const NT a, const NT b) { return a + b; });
653
654 return result_sum + result_integral;
655 }
656 } // namespace LoopIntegrals
657} // namespace DiFfRG
DiFfRG::LoopIntegrals::spatial_angle_integrate_and_sum
NT spatial_angle_integrate_and_sum(const FUN &fun, const QGauss< 1 > &x_quadrature, const double x_extent, const double k, const QGauss< 1 > &cos_quadrature, const int m_order, const QGauss< 1 > &m_quadrature, const double m_extent, const double T)
Performs the integral.
Definition loop_integrals.hh:529
DiFfRG::LoopIntegrals::two_angle_integrate
NT two_angle_integrate(const FUN &fun, const QGauss< 1 > &x_quadrature, const double x_extent, const double k, const QGauss< 1 > &cos_quadrature)
Performs the integral.
Definition loop_integrals.hh:145
DiFfRG::LoopIntegrals::spatial_integrate_and_integrate
NT spatial_integrate_and_integrate(const FUN &fun, const QGauss< 1 > &x_quadrature, const double x_extent, const double k, const QGauss< 1 > &m_quadrature, const double m_extent)
Performs the integral.
Definition loop_integrals.hh:296
DiFfRG::LoopIntegrals::spatial_angle_integrate_and_integrate
NT spatial_angle_integrate_and_integrate(const FUN &fun, const QGauss< 1 > &x_quadrature, const double x_extent, const double k, const QGauss< 1 > &cos_quadrature, const QGauss< 1 > &m_quadrature, const double m_extent)
Performs the integral.
Definition loop_integrals.hh:360
DiFfRG::LoopIntegrals::angle_integrate
NT angle_integrate(const FUN &fun, const QGauss< 1 > &x_quadrature, const double x_extent, const double k, const QGauss< 1 > &cos_quadrature)
Performs the integral.
Definition loop_integrals.hh:84
DiFfRG::LoopIntegrals::three_angle_integrate
NT three_angle_integrate(const FUN &fun, const QGauss< 1 > &x_quadrature, const double x_extent, const double k, const QGauss< 1 > &cos_quadrature)
Performs the integral.
Definition loop_integrals.hh:215
DiFfRG::LoopIntegrals::integrate
NT integrate(const FUN &fun, const QGauss< 1 > &x_quadrature, const double x_extent, const double k)
Performs the integral.
Definition loop_integrals.hh:35
DiFfRG::LoopIntegrals::sum
NT sum(const FUN &fun, const int m_order, const QGauss< 1 > &m_quadrature, const double m_extent, const double T)
Definition loop_integrals.hh:619
DiFfRG::LoopIntegrals::spatial_integrate_and_sum
NT spatial_integrate_and_sum(const FUN &fun, const QGauss< 1 > &x_quadrature, const double x_extent, const double k, const int m_order, const QGauss< 1 > &m_quadrature, const double m_extent, const double T)
Performs the integral.
Definition loop_integrals.hh:436
DiFfRG
Definition complex_math.hh:14
DiFfRG::powr
constexpr __forceinline__ __host__ __device__ NumberType powr(const NumberType x)
A compile-time evaluatable power function for whole number exponents.
Definition math.hh:45
DiFfRG::is_close
bool __forceinline__ __host__ __device__ is_close(T1 a, T2 b, T3 eps_)
Function to evaluate whether two floats are equal to numerical precision. Tests for both relative and...
Definition math.hh:160
DiFfRG::S_d
constexpr __forceinline__ __host__ __device__ double S_d(NT d)
Surface of a d-dimensional sphere.
Definition math.hh:91