Block-Structured AMR Software Framework
AMReX_Particle_mod_K.H
Go to the documentation of this file.
1 #ifndef AMREX_PARTICLE_MOD_K_H_
2 #define AMREX_PARTICLE_MOD_K_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_FArrayBox.H>
6 #include <AMReX_REAL.H>
7 #include <cmath>
8 
9 namespace amrex {
10 
11 template <typename P>
12 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
13 void amrex_deposit_cic (P const& p, int nc, amrex::Array4<amrex::Real> const& rho,
14  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
15  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi)
16 {
17 // GCC does not like rdata(comp)
18 #if defined(__GNUC__) && !defined(__clang__)
19 #pragma GCC diagnostic push
20 #pragma GCC diagnostic ignored "-Warray-bounds"
21 #endif
22 
23 #if (AMREX_SPACEDIM == 1)
24  amrex::Real lx = (p.pos(0) - plo[0]) * dxi[0] + Real(0.5);
25 
26  int i = static_cast<int>(amrex::Math::floor(lx));
27 
28  amrex::Real xint = lx - static_cast<Real>(i);
29 
30  amrex::Real sx[2] = {Real(1.0) - xint, xint};
31 
32  for (int ii = 0; ii <= 1; ++ii) {
33  amrex::Gpu::Atomic::AddNoRet(&rho(i+ii-1, 0, 0, 0), static_cast<Real>(sx[ii]*p.rdata(0)));
34  }
35 
36  for (int comp=1; comp < nc; ++comp) {
37  for (int ii = 0; ii <= 1; ++ii) {
38  amrex::Gpu::Atomic::AddNoRet(&rho(i+ii-1, 0, 0, comp),
39  static_cast<Real>(sx[ii]*p.rdata(0)*p.rdata(comp)));
40  }
41  }
42 #elif (AMREX_SPACEDIM == 2)
43  amrex::Real lx = (p.pos(0) - plo[0]) * dxi[0] + Real(0.5);
44  amrex::Real ly = (p.pos(1) - plo[1]) * dxi[1] + Real(0.5);
45 
46  int i = static_cast<int>(amrex::Math::floor(lx));
47  int j = static_cast<int>(amrex::Math::floor(ly));
48 
49  amrex::Real xint = lx - static_cast<Real>(i);
50  amrex::Real yint = ly - static_cast<Real>(j);
51 
52  amrex::Real sx[2] = {Real(1.0) - xint, xint};
53  amrex::Real sy[2] = {Real(1.0) - yint, yint};
54 
55  for (int jj = 0; jj <= 1; ++jj) {
56  for (int ii = 0; ii <= 1; ++ii) {
57  amrex::Gpu::Atomic::AddNoRet(&rho(i+ii-1, j+jj-1, 0, 0),
58  static_cast<Real>(sx[ii]*sy[jj]*p.rdata(0)));
59  }
60  }
61 
62  for (int comp=1; comp < nc; ++comp) {
63  for (int jj = 0; jj <= 1; ++jj) {
64  for (int ii = 0; ii <= 1; ++ii) {
65  amrex::Gpu::Atomic::AddNoRet(&rho(i+ii-1, j+jj-1, 0, comp),
66  static_cast<Real>(sx[ii]*sy[jj]*p.rdata(0)*p.rdata(comp)));
67  }
68  }
69  }
70 
71 #elif (AMREX_SPACEDIM == 3)
72  amrex::Real lx = (p.pos(0) - plo[0]) * dxi[0] + Real(0.5);
73  amrex::Real ly = (p.pos(1) - plo[1]) * dxi[1] + Real(0.5);
74  amrex::Real lz = (p.pos(2) - plo[2]) * dxi[2] + Real(0.5);
75 
76  int i = static_cast<int>(amrex::Math::floor(lx));
77  int j = static_cast<int>(amrex::Math::floor(ly));
78  int k = static_cast<int>(amrex::Math::floor(lz));
79 
80  amrex::Real xint = lx - static_cast<Real>(i);
81  amrex::Real yint = ly - static_cast<Real>(j);
82  amrex::Real zint = lz - static_cast<Real>(k);
83 
84  amrex::Real sx[] = {Real(1.0) - xint, xint};
85  amrex::Real sy[] = {Real(1.0) - yint, yint};
86  amrex::Real sz[] = {Real(1.0) - zint, zint};
87 
88  for (int kk = 0; kk <= 1; ++kk) {
89  for (int jj = 0; jj <= 1; ++jj) {
90  for (int ii = 0; ii <= 1; ++ii) {
91  amrex::Gpu::Atomic::AddNoRet(&rho(i+ii-1, j+jj-1, k+kk-1, 0),
92  static_cast<Real>(sx[ii]*sy[jj]*sz[kk]*p.rdata(0)));
93  }
94  }
95  }
96 
97  for (int comp=1; comp < nc; ++comp) {
98  for (int kk = 0; kk <= 1; ++kk) {
99  for (int jj = 0; jj <= 1; ++jj) {
100  for (int ii = 0; ii <= 1; ++ii) {
101  amrex::Gpu::Atomic::AddNoRet(&rho(i+ii-1, j+jj-1, k+kk-1, comp),
102  static_cast<Real>(sx[ii]*sy[jj]*sz[kk]*p.rdata(0)*p.rdata(comp)));
103  }
104  }
105  }
106  }
107 #else
108  amrex::Abort("Not implemented.");
109 #endif
110 
111 #if defined(__GNUC__) && !defined(__clang__)
112 #pragma GCC diagnostic pop
113 #endif
114 }
115 
116 template <typename P>
117 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
118 void amrex_deposit_particle_dx_cic (P const& p, int nc, amrex::Array4<amrex::Real> const& rho,
119  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
120  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
121  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& pdxi)
122 {
123 // GCC does not like rdata(comp)
124 #if defined(__GNUC__) && !defined(__clang__)
125 #pragma GCC diagnostic push
126 #pragma GCC diagnostic ignored "-Warray-bounds"
127 #endif
128 
129 #if (AMREX_SPACEDIM == 1)
130  amrex::Real factor = (pdxi[0]/dxi[0]);
131 
132  amrex::Real lx = (p.pos(0) - plo[0] - Real(0.5)/pdxi[0]) * dxi[0];
133 
134  amrex::Real hx = (p.pos(0) - plo[0] + Real(0.5)/pdxi[0]) * dxi[0];
135 
136  int lo_x = static_cast<int>(amrex::Math::floor(lx));
137 
138  int hi_x = static_cast<int>(amrex::Math::floor(hx));
139 
140  for (int i = lo_x; i <= hi_x; ++i) {
141  if (i < rho.begin.x || i >= rho.end.x) { continue; }
142  amrex::Real wx = amrex::min(hx - static_cast<Real>(i), amrex::Real(1.0)) - amrex::max(lx - static_cast<Real>(i), amrex::Real(0.0));
143  amrex::Real weight = wx*factor;
144  amrex::Gpu::Atomic::AddNoRet(&rho(i, 0, 0, 0), static_cast<Real>(weight*p.rdata(0)));
145  }
146 
147  for (int comp = 1; comp < nc; ++comp)
148  {
149  for (int i = lo_x; i <= hi_x; ++i) {
150  if (i < rho.begin.x || i >= rho.end.x) { continue; }
151  amrex::Real wx = amrex::min(hx - static_cast<Real>(i), amrex::Real(1.0)) - amrex::max(lx - static_cast<Real>(i), amrex::Real(0.0));
152  amrex::Real weight = wx*factor;
153  amrex::Gpu::Atomic::AddNoRet(&rho(i, 0, 0, comp), static_cast<Real>(weight*p.rdata(0)*p.rdata(comp)));
154  }
155  }
156 
157 #elif (AMREX_SPACEDIM == 2)
158  amrex::Real factor = (pdxi[0]/dxi[0])*(pdxi[1]/dxi[1]);
159 
160  amrex::Real lx = (p.pos(0) - plo[0] - Real(0.5)/pdxi[0]) * dxi[0];
161  amrex::Real ly = (p.pos(1) - plo[1] - Real(0.5)/pdxi[1]) * dxi[1];
162 
163  amrex::Real hx = (p.pos(0) - plo[0] + Real(0.5)/pdxi[0]) * dxi[0];
164  amrex::Real hy = (p.pos(1) - plo[1] + Real(0.5)/pdxi[1]) * dxi[1];
165 
166  int lo_x = static_cast<int>(amrex::Math::floor(lx));
167  int lo_y = static_cast<int>(amrex::Math::floor(ly));
168 
169  int hi_x = static_cast<int>(amrex::Math::floor(hx));
170  int hi_y = static_cast<int>(amrex::Math::floor(hy));
171 
172  for (int j = lo_y; j <= hi_y; ++j) {
173  if (j < rho.begin.y || j >= rho.end.y) { continue; }
174  amrex::Real wy = amrex::min(hy - static_cast<Real>(j), amrex::Real(1.0)) - amrex::max(ly - static_cast<Real>(j), amrex::Real(0.0));
175  for (int i = lo_x; i <= hi_x; ++i) {
176  if (i < rho.begin.x || i >= rho.end.x) { continue; }
177  amrex::Real wx = amrex::min(hx - static_cast<Real>(i), amrex::Real(1.0)) - amrex::max(lx - static_cast<Real>(i), amrex::Real(0.0));
178  amrex::Real weight = wx*wy*factor;
179  amrex::Gpu::Atomic::AddNoRet(&rho(i, j, 0, 0), static_cast<Real>(weight*p.rdata(0)));
180  }
181  }
182 
183  for (int comp = 1; comp < nc; ++comp) {
184  for (int j = lo_y; j <= hi_y; ++j) {
185  if (j < rho.begin.y || j >= rho.end.y) { continue; }
186  amrex::Real wy = amrex::min(hy - static_cast<Real>(j), amrex::Real(1.0)) - amrex::max(ly - static_cast<Real>(j), amrex::Real(0.0));
187  for (int i = lo_x; i <= hi_x; ++i) {
188  if (i < rho.begin.x || i >= rho.end.x) { continue; }
189  amrex::Real wx = amrex::min(hx - static_cast<Real>(i), amrex::Real(1.0)) - amrex::max(lx - static_cast<Real>(i), amrex::Real(0.0));
190  amrex::Real weight = wx*wy*factor;
191  amrex::Gpu::Atomic::AddNoRet(&rho(i, j, 0, comp), static_cast<Real>(weight*p.rdata(0)*p.rdata(comp)));
192  }
193  }
194  }
195 
196 #elif (AMREX_SPACEDIM == 3)
197  amrex::Real factor = (pdxi[0]/dxi[0])*(pdxi[1]/dxi[1])*(pdxi[2]/dxi[2]);
198 
199  amrex::Real lx = (p.pos(0) - plo[0] - Real(0.5)/pdxi[0]) * dxi[0];
200  amrex::Real ly = (p.pos(1) - plo[1] - Real(0.5)/pdxi[1]) * dxi[1];
201  amrex::Real lz = (p.pos(2) - plo[2] - Real(0.5)/pdxi[2]) * dxi[2];
202 
203  amrex::Real hx = (p.pos(0) - plo[0] + Real(0.5)/pdxi[0]) * dxi[0];
204  amrex::Real hy = (p.pos(1) - plo[1] + Real(0.5)/pdxi[1]) * dxi[1];
205  amrex::Real hz = (p.pos(2) - plo[2] + Real(0.5)/pdxi[2]) * dxi[2];
206 
207  int lo_x = static_cast<int>(amrex::Math::floor(lx));
208  int lo_y = static_cast<int>(amrex::Math::floor(ly));
209  int lo_z = static_cast<int>(amrex::Math::floor(lz));
210 
211  int hi_x = static_cast<int>(amrex::Math::floor(hx));
212  int hi_y = static_cast<int>(amrex::Math::floor(hy));
213  int hi_z = static_cast<int>(amrex::Math::floor(hz));
214 
215  for (int k = lo_z; k <= hi_z; ++k) {
216  if (k < rho.begin.z || k >= rho.end.z) { continue; }
217  amrex::Real wz = amrex::min(hz - static_cast<Real>(k), amrex::Real(1.0)) - amrex::max(lz - static_cast<Real>(k), amrex::Real(0.0));
218  for (int j = lo_y; j <= hi_y; ++j) {
219  if (j < rho.begin.y || j >= rho.end.y) { continue; }
220  amrex::Real wy = amrex::min(hy - static_cast<Real>(j), amrex::Real(1.0)) - amrex::max(ly - static_cast<Real>(j), amrex::Real(0.0));
221  for (int i = lo_x; i <= hi_x; ++i) {
222  if (i < rho.begin.x || i >= rho.end.x) { continue; }
223  amrex::Real wx = amrex::min(hx - static_cast<Real>(i), amrex::Real(1.0)) - amrex::max(lx - static_cast<Real>(i), amrex::Real(0.0));
224  amrex::Real weight = wx*wy*wz*factor;
225  amrex::Gpu::Atomic::AddNoRet(&rho(i, j, k, 0), static_cast<Real>(weight*p.rdata(0)));
226  }
227  }
228  }
229 
230  for (int comp = 1; comp < nc; ++comp) {
231  for (int k = lo_z; k <= hi_z; ++k) {
232  if (k < rho.begin.z || k >= rho.end.z) { continue; }
233  amrex::Real wz = amrex::min(hz - static_cast<Real>(k), amrex::Real(1.0)) - amrex::max(lz - static_cast<Real>(k), amrex::Real(0.0));
234  for (int j = lo_y; j <= hi_y; ++j) {
235  if (j < rho.begin.y || j >= rho.end.y) { continue; }
236  amrex::Real wy = amrex::min(hy - static_cast<Real>(j), amrex::Real(1.0)) - amrex::max(ly - static_cast<Real>(j), amrex::Real(0.0));
237  for (int i = lo_x; i <= hi_x; ++i) {
238  if (i < rho.begin.x || i >= rho.end.x) { continue; }
239  amrex::Real wx = amrex::min(hx - static_cast<Real>(i), amrex::Real(1.0)) - amrex::max(lx - static_cast<Real>(i), amrex::Real(0.0));
240  amrex::Real weight = wx*wy*wz*factor;
241  amrex::Gpu::Atomic::AddNoRet(&rho(i, j, k, comp), static_cast<Real>(weight*p.rdata(0)*p.rdata(comp)));
242  }
243  }
244  }
245  }
246 #else
247  amrex::Abort("Not implemented.")
248 #endif
249 
250 #if defined(__GNUC__) && !defined(__clang__)
251 #pragma GCC diagnostic pop
252 #endif
253 }
254 
255 }
256 
257 #endif
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
AMReX_FArrayBox.H
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
AMReX_REAL.H
amrex::Gpu::Atomic::AddNoRet
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AddNoRet(T *sum, T value) noexcept
Definition: AMReX_GpuAtomic.H:281
amrex::openbc::P
static constexpr int P
Definition: AMReX_OpenBC.H:14
amrex
Definition: AMReX_Amr.cpp:49
amrex::amrex_deposit_cic
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void amrex_deposit_cic(P const &p, int nc, amrex::Array4< amrex::Real > const &rho, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &plo, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &dxi)
Definition: AMReX_Particle_mod_K.H:13
amrex::max
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
amrex::min
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:21
amrex::amrex_deposit_particle_dx_cic
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void amrex_deposit_particle_dx_cic(P const &p, int nc, amrex::Array4< amrex::Real > const &rho, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &plo, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &dxi, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &pdxi)
Definition: AMReX_Particle_mod_K.H:118
amrex::Abort
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition: AMReX.cpp:225
amrex::Array4
Definition: AMReX_Array4.H:61
amrex::Array4::end
Dim3 end
Definition: AMReX_Array4.H:67
amrex::Dim3::x
int x
Definition: AMReX_Dim3.H:12
amrex::Dim3::z
int z
Definition: AMReX_Dim3.H:12
amrex::Dim3::y
int y
Definition: AMReX_Dim3.H:12
amrex::GpuArray
Definition: AMReX_Array.H:34