ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
math.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2024 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22#pragma once
23
25#include <utils/math/sqr.hpp>
26
27#ifndef __CUDACC__
28#include <cmath>
29#endif
30#include <numbers>
31#include <stdexcept>
32#include <string>
33
34namespace math {
35
36/** @brief Return the absolute value of x. */
37inline DEVICE_QUALIFIER auto abs(double x) { return fabs(x); }
38
39/** @brief Return the absolute value of x. */
40inline DEVICE_QUALIFIER auto abs(float x) { return fabsf(x); }
41
42/**
43 * @brief Calculate the function @f$ \mathrm{sinc}(x) = \sin(\pi x)/(\pi x) @f$.
44 *
45 * (same convention as in @cite hockney88a). In order to avoid divisions
46 * by 0, arguments whose modulus is smaller than @f$ \epsilon = 0.1 @f$
47 * are evaluated by an 8th order Taylor expansion.
48 * Note that the difference between sinc(x) and this expansion
49 * is smaller than 0.235e-12, if x is smaller than @f$ \epsilon @f$.
50 * (The next term in the expansion is the 10th order contribution, i.e.
51 * @f$ \pi^{10} x^{10}/39916800 \approx 0.2346 \cdot x^{12} @f$).
52 */
53template <typename T> DEVICE_QUALIFIER auto sinc(T x) {
54 auto constexpr epsilon = T(0.1);
55#if not defined(__CUDACC__)
56 using std::sin;
57#endif
58
59 auto const pix = std::numbers::pi_v<T> * x;
60
61 if (::math::abs(x) > epsilon)
62 return sin(pix) / pix;
63
64 auto constexpr factorial = [](int n) consteval {
65 int acc{1}, c{1};
66 while (c < n) {
67 acc *= ++c;
68 }
69 return acc;
70 };
71
72 /* Coefficients of the Taylor expansion of sinc */
73 auto constexpr c0 = T(+1) / T(factorial(1));
74 auto constexpr c2 = T(-1) / T(factorial(3));
75 auto constexpr c4 = T(+1) / T(factorial(5));
76 auto constexpr c6 = T(-1) / T(factorial(7));
77 auto constexpr c8 = T(+1) / T(factorial(9));
78
79 auto const pix2 = pix * pix;
80 return c0 + pix2 * (c2 + pix2 * (c4 + pix2 * (c6 + pix2 * c8)));
81}
82
83/**
84 * @brief One of the aliasing sums used to compute k-space errors.
85 * Fortunately the one which is most important (because it converges
86 * most slowly, since it is not damped exponentially) can be calculated
87 * analytically. The result (which depends on the order of the spline
88 * interpolation) can be written as an even trigonometric polynomial.
89 * The results are tabulated here (the employed formula is eq. 7-66
90 * p. 233 in @cite hockney88a).
91 */
92template <int cao>
93DEVICE_QUALIFIER auto analytic_cotangent_sum(int n, double mesh_i) {
94 static_assert(cao >= 1 and cao <= 7);
95#if not defined(__CUDACC__)
96 using std::cos;
97#endif
98 auto const theta = static_cast<double>(n) * mesh_i * std::numbers::pi;
99 auto const c = Utils::sqr(cos(theta));
100
101 if constexpr (cao == 1) {
102 return 1.;
103 }
104 if constexpr (cao == 2) {
105 return (1. + c * 2.) / 3.;
106 }
107 if constexpr (cao == 3) {
108 return (2. + c * (11. + c * 2.)) / 15.;
109 }
110 if constexpr (cao == 4) {
111 return (17. + c * (180. + c * (114. + c * 4.))) / 315.;
112 }
113 if constexpr (cao == 5) {
114 return (62. + c * (1072. + c * (1452. + c * (247. + c * 2.)))) / 2835.;
115 }
116 if constexpr (cao == 6) {
117 return (1382. +
118 c * (35396. + c * (83021. + c * (34096. + c * (2026. + c * 4.))))) /
119 155925.;
120 }
121 return (21844. +
122 c * (776661. +
123 c * (2801040. +
124 c * (2123860. + c * (349500. + c * (8166. + c * 4.)))))) /
125 6081075.;
126}
127
129 decltype(&analytic_cotangent_sum<1>) ptr = nullptr;
130 if (cao == 1) {
131 ptr = &analytic_cotangent_sum<1>;
132 } else if (cao == 2) {
133 ptr = &analytic_cotangent_sum<2>;
134 } else if (cao == 3) {
135 ptr = &analytic_cotangent_sum<3>;
136 } else if (cao == 4) {
137 ptr = &analytic_cotangent_sum<4>;
138 } else if (cao == 5) {
139 ptr = &analytic_cotangent_sum<5>;
140 } else if (cao == 6) {
141 ptr = &analytic_cotangent_sum<6>;
142 } else if (cao == 7) {
143 ptr = &analytic_cotangent_sum<7>;
144 }
145 if (ptr == nullptr) {
146 throw std::logic_error("Invalid value cao=" + std::to_string(cao));
147 }
148 return ptr;
149}
150
151} // namespace math
#define DEVICE_QUALIFIER
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
Definition math.hpp:34
DEVICE_QUALIFIER auto abs(double x)
Return the absolute value of x.
Definition math.hpp:37
DEVICE_QUALIFIER auto sinc(T x)
Calculate the function .
Definition math.hpp:53
auto get_analytic_cotangent_sum_kernel(int cao)
Definition math.hpp:128
DEVICE_QUALIFIER auto analytic_cotangent_sum(int n, double mesh_i)
One of the aliasing sums used to compute k-space errors.
Definition math.hpp:93