ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
p3m_gpu_error_cuda.cu
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19/** @file
20 *
21 * P3M electrostatics on GPU.
22 *
23 * The corresponding header file is p3m_gpu_error.hpp.
24 */
25
26#include "p3m/math.hpp"
27#include "p3m_gpu_error.hpp"
28
29#include "cuda/utils.cuh"
30
32#include <utils/math/sqr.hpp>
33
34#include <thrust/device_vector.h>
35#include <thrust/reduce.h>
36
37#include <cuda.h>
38
39#include <algorithm>
40#include <cstddef>
41#include <numbers>
42
43#if defined(OMPI_MPI_H) || defined(_MPI_H)
44#error CU-file includes mpi.h! This should not happen!
45#endif
46
47template <int cao>
48__global__ void p3m_k_space_error_gpu_kernel_ik(int3 mesh, double3 meshi,
49 double alpha_L, double *he_q) {
50 using Utils::int_pow;
51 using Utils::sqr;
52
53 const int nx =
54 -mesh.x / 2 + static_cast<int>(blockDim.x * blockIdx.x + threadIdx.x);
55 const int ny =
56 -mesh.y / 2 + static_cast<int>(blockDim.y * blockIdx.y + threadIdx.y);
57 const int nz =
58 -mesh.z / 2 + static_cast<int>(blockDim.z * blockIdx.z + threadIdx.z);
59
60 if ((nx >= mesh.x / 2) || (ny >= mesh.y / 2) || (nz >= mesh.z / 2))
61 return;
62
63 const int lind = (nx + mesh.x / 2) * mesh.y * mesh.z +
64 (ny + mesh.y / 2) * mesh.z + (nz + mesh.z / 2);
65
66 if ((nx != 0) || (ny != 0) || (nz != 0)) {
67 const double alpha_L_i = 1. / alpha_L;
68 const double n2 = sqr(nx) + sqr(ny) + sqr(nz);
69 const double cs = math::analytic_cotangent_sum<cao>(nz, meshi.z) *
70 math::analytic_cotangent_sum<cao>(nx, meshi.x) *
71 math::analytic_cotangent_sum<cao>(ny, meshi.y);
72 const double ex = exp(-sqr(std::numbers::pi * alpha_L_i) * n2);
73 const double ex2 = sqr(ex);
74 const double U2 =
75 int_pow<2 * cao>(math::sinc(meshi.x * nx) * math::sinc(meshi.y * ny) *
76 math::sinc(meshi.z * nz));
77 auto const alias1 = ex2 / n2;
78 auto const d = alias1 - sqr(U2 * ex / cs) / n2;
79
80 if (d > 0. and (d / alias1 > ROUND_ERROR_PREC))
81 he_q[lind] = d;
82 } else {
83 he_q[lind] = 0.;
84 }
85}
86
87double p3m_k_space_error_gpu(double prefactor, const int *mesh, int cao,
88 int npart, double sum_q2, double alpha_L,
89 const double *box) {
90 static thrust::device_vector<double> he_q;
91 he_q.resize(static_cast<unsigned>(mesh[0] * mesh[1] * mesh[2]));
92
93 dim3 grid(std::max<unsigned>(1u, static_cast<unsigned>(mesh[0] / 8 + 1)),
94 std::max<unsigned>(1u, static_cast<unsigned>(mesh[1] / 8 + 1)),
95 std::max<unsigned>(1u, static_cast<unsigned>(mesh[2] / 8 + 1)));
96
97 dim3 block(8, 8, 8);
98
99 int3 mesh3;
100 mesh3.x = mesh[0];
101 mesh3.y = mesh[1];
102 mesh3.z = mesh[2];
103
104 double3 meshi;
105 meshi.x = 1. / mesh[0];
106 meshi.y = 1. / mesh[1];
107 meshi.z = 1. / mesh[2];
108
109 switch (cao) {
110 case 1:
111 KERNELCALL(p3m_k_space_error_gpu_kernel_ik<1>, grid, block, mesh3, meshi,
112 alpha_L, thrust::raw_pointer_cast(he_q.data()));
113 break;
114 case 2:
115 KERNELCALL(p3m_k_space_error_gpu_kernel_ik<2>, grid, block, mesh3, meshi,
116 alpha_L, thrust::raw_pointer_cast(he_q.data()));
117 break;
118 case 3:
119 KERNELCALL(p3m_k_space_error_gpu_kernel_ik<3>, grid, block, mesh3, meshi,
120 alpha_L, thrust::raw_pointer_cast(he_q.data()));
121 break;
122 case 4:
123 KERNELCALL(p3m_k_space_error_gpu_kernel_ik<4>, grid, block, mesh3, meshi,
124 alpha_L, thrust::raw_pointer_cast(he_q.data()));
125 break;
126 case 5:
127 KERNELCALL(p3m_k_space_error_gpu_kernel_ik<5>, grid, block, mesh3, meshi,
128 alpha_L, thrust::raw_pointer_cast(he_q.data()));
129 break;
130 case 6:
131 KERNELCALL(p3m_k_space_error_gpu_kernel_ik<6>, grid, block, mesh3, meshi,
132 alpha_L, thrust::raw_pointer_cast(he_q.data()));
133 break;
134 case 7:
135 KERNELCALL(p3m_k_space_error_gpu_kernel_ik<7>, grid, block, mesh3, meshi,
136 alpha_L, thrust::raw_pointer_cast(he_q.data()));
137 break;
138 }
139
140 auto const he_q_final = thrust::reduce(he_q.begin(), he_q.end());
141
142 return 2. * prefactor * sum_q2 * sqrt(he_q_final / npart) / (box[1] * box[2]);
143}
#define ROUND_ERROR_PREC
Precision for capture of round off errors.
Definition config.hpp:66
static double * block(double *p, std::size_t index, std::size_t size)
Definition elc.cpp:172
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
DEVICE_QUALIFIER constexpr T int_pow(T x)
Calculate integer powers.
Definition int_pow.hpp:61
DEVICE_QUALIFIER auto sinc(T x)
Calculate the function .
Definition math.hpp:53
P3M electrostatics on GPU.
double p3m_k_space_error_gpu(double prefactor, const int *mesh, int cao, int npart, double sum_q2, double alpha_L, const double *box)
__global__ void p3m_k_space_error_gpu_kernel_ik(int3 mesh, double3 meshi, double alpha_L, double *he_q)
#define KERNELCALL(_function, _grid, _block,...)
Definition utils.cuh:79