ARTEMIS
WarpXComm_K.H
Go to the documentation of this file.
1 /* Copyright 2019 Axel Huebl, Weiqun Zhang
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef WARPX_COMM_K_H_
8 #define WARPX_COMM_K_H_
9 
10 #include <AMReX.H>
11 #include <AMReX_FArrayBox.H>
12 
13 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
14 void warpx_interp (int j, int k, int l,
15  amrex::Array4<amrex::Real > const& arr_aux,
16  amrex::Array4<amrex::Real const> const& arr_fine,
17  amrex::Array4<amrex::Real const> const& arr_coarse,
18  const amrex::IntVect& arr_stag,
19  const amrex::IntVect& rr)
20 {
21  using namespace amrex;
22 
23  // Pad arr_coarse with zeros beyond ghost cells for out-of-bound accesses
24  const auto arr_coarse_zeropad = [arr_coarse] (const int jj, const int kk, const int ll) noexcept
25  {
26  return arr_coarse.contains(jj,kk,ll) ? arr_coarse(jj,kk,ll) : 0.0_rt;
27  };
28 
29  // NOTE Indices (j,k,l) in the following refer to (z,-,-) in 1D, (x,z,-) in 2D, and (x,y,z) in 3D
30 
31  // Refinement ratio
32  const int rj = rr[0];
33  const int rk = (AMREX_SPACEDIM == 1) ? 1 : rr[1];
34  const int rl = (AMREX_SPACEDIM <= 2) ? 1 : rr[2];
35 
36  // Staggering (0: cell-centered; 1: nodal)
37  const int sj = arr_stag[0];
38  const int sk = (AMREX_SPACEDIM == 1) ? 0 : arr_stag[1];
39  const int sl = (AMREX_SPACEDIM <= 2) ? 0 : arr_stag[2];
40 
41  // Number of points used for interpolation from coarse grid to fine grid
42  const int nj = 2;
43  const int nk = 2;
44  const int nl = 2;
45 
46  const int jc = (sj == 0) ? amrex::coarsen(j - rj/2, rj) : amrex::coarsen(j, rj);
47  const int kc = (sk == 0) ? amrex::coarsen(k - rk/2, rk) : amrex::coarsen(k, rk);
48  const int lc = (sl == 0) ? amrex::coarsen(l - rl/2, rl) : amrex::coarsen(l, rl);
49 
50  amrex::Real wj;
51  amrex::Real wk;
52  amrex::Real wl;
53 
54  // Interpolate from coarse grid to fine grid using 2 points
55  // with weights depending on the distance, for both nodal and cell-centered grids
56  amrex::Real hj = (sj == 0) ? 0.5_rt : 0._rt;
57  amrex::Real hk = (sk == 0) ? 0.5_rt : 0._rt;
58  amrex::Real hl = (sl == 0) ? 0.5_rt : 0._rt;
59 
60  amrex::Real res = 0.0_rt;
61 
62  for (int jj = 0; jj < nj; jj++) {
63  for (int kk = 0; kk < nk; kk++) {
64  for (int ll = 0; ll < nl; ll++) {
65  wj = (rj - amrex::Math::abs(j + hj - (jc + jj + hj) * rj)) / static_cast<amrex::Real>(rj);
66  wk = (rk - amrex::Math::abs(k + hk - (kc + kk + hk) * rk)) / static_cast<amrex::Real>(rk);
67  wl = (rl - amrex::Math::abs(l + hl - (lc + ll + hl) * rl)) / static_cast<amrex::Real>(rl);
68  res += wj * wk * wl * arr_coarse_zeropad(jc+jj,kc+kk,lc+ll);
69  }
70  }
71  }
72  arr_aux(j,k,l) = arr_fine(j,k,l) + res;
73 }
74 
75 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
76 void warpx_interp_nd_bfield_x (int j, int k, int l,
77  amrex::Array4<amrex::Real> const& Bxa,
78  amrex::Array4<amrex::Real const> const& Bxf,
79  amrex::Array4<amrex::Real const> const& Bxc,
80  amrex::Array4<amrex::Real const> const& Bxg)
81 {
82  using namespace amrex;
83 
84  // Pad Bxf with zeros beyond ghost cells for out-of-bound accesses
85  const auto Bxf_zeropad = [Bxf] (const int jj, const int kk, const int ll) noexcept
86  {
87  return Bxf.contains(jj,kk,ll) ? Bxf(jj,kk,ll) : 0.0_rt;
88  };
89 
90  int jg = amrex::coarsen(j,2);
91  Real wx = (j == jg*2) ? 0.0_rt : 0.5_rt;
92  Real owx = 1.0_rt-wx;
93 
94  int kg = amrex::coarsen(k,2);
95  Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
96  Real owy = 1.0_rt-wy;
97 
98 #if defined(WARPX_DIM_1D_Z)
99 
100  // interp from coarse nodal to fine nodal
101  Real bg = owx * Bxg(jg ,0,0)
102  + wx * Bxg(jg+1,0,0);
103 
104  // interp from coarse staggered to fine nodal
105  Real bc = owx * Bxc(jg ,0,0)
106  + wx * Bxc(jg+1,0,0);
107 
108  // interp from fine staggered to fine nodal
109  Real bf = 0.5_rt*(Bxf_zeropad(j-1,0,0) + Bxf_zeropad(j,0,0));
110  amrex::ignore_unused(owy);
111 
112 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
113 
114  // interp from coarse nodal to fine nodal
115  Real bg = owx * owy * Bxg(jg ,kg ,0)
116  + owx * wy * Bxg(jg ,kg+1,0)
117  + wx * owy * Bxg(jg+1,kg ,0)
118  + wx * wy * Bxg(jg+1,kg+1,0);
119 
120  // interp from coarse staggered to fine nodal
121  wy = 0.5_rt-wy; owy = 1.0_rt-wy;
122  Real bc = owx * owy * Bxc(jg ,kg ,0)
123  + owx * wy * Bxc(jg ,kg-1,0)
124  + wx * owy * Bxc(jg+1,kg ,0)
125  + wx * wy * Bxc(jg+1,kg-1,0);
126 
127  // interp from fine staggered to fine nodal
128  Real bf = 0.5_rt*(Bxf_zeropad(j,k-1,0) + Bxf_zeropad(j,k,0));
129 
130 #else
131 
132  int lg = amrex::coarsen(l,2);
133  Real wz = (l == lg*2) ? 0.0_rt : 0.5_rt;
134  Real owz = 1.0_rt-wz;
135 
136  // interp from coarse nodal to fine nodal
137  Real bg = owx * owy * owz * Bxg(jg ,kg ,lg )
138  + wx * owy * owz * Bxg(jg+1,kg ,lg )
139  + owx * wy * owz * Bxg(jg ,kg+1,lg )
140  + wx * wy * owz * Bxg(jg+1,kg+1,lg )
141  + owx * owy * wz * Bxg(jg ,kg ,lg+1)
142  + wx * owy * wz * Bxg(jg+1,kg ,lg+1)
143  + owx * wy * wz * Bxg(jg ,kg+1,lg+1)
144  + wx * wy * wz * Bxg(jg+1,kg+1,lg+1);
145 
146  // interp from coarse staggered to fine nodal
147  wy = 0.5_rt-wy; owy = 1.0_rt-wy;
148  wz = 0.5_rt-wz; owz = 1.0_rt-wz;
149  Real bc = owx * owy * owz * Bxc(jg ,kg ,lg )
150  + wx * owy * owz * Bxc(jg+1,kg ,lg )
151  + owx * wy * owz * Bxc(jg ,kg-1,lg )
152  + wx * wy * owz * Bxc(jg+1,kg-1,lg )
153  + owx * owy * wz * Bxc(jg ,kg ,lg-1)
154  + wx * owy * wz * Bxc(jg+1,kg ,lg-1)
155  + owx * wy * wz * Bxc(jg ,kg-1,lg-1)
156  + wx * wy * wz * Bxc(jg+1,kg-1,lg-1);
157 
158  // interp from fine stagged to fine nodal
159  Real bf = 0.25_rt*(Bxf_zeropad(j,k-1,l-1) + Bxf_zeropad(j,k,l-1) + Bxf_zeropad(j,k-1,l) + Bxf_zeropad(j,k,l));
160 #endif
161 
162  Bxa(j,k,l) = bg + (bf-bc);
163 }
164 
165 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
166 void warpx_interp_nd_bfield_y (int j, int k, int l,
167  amrex::Array4<amrex::Real> const& Bya,
168  amrex::Array4<amrex::Real const> const& Byf,
169  amrex::Array4<amrex::Real const> const& Byc,
170  amrex::Array4<amrex::Real const> const& Byg)
171 {
172  using namespace amrex;
173 
174  // Pad Byf with zeros beyond ghost cells for out-of-bound accesses
175  const auto Byf_zeropad = [Byf] (const int jj, const int kk, const int ll) noexcept
176  {
177  return Byf.contains(jj,kk,ll) ? Byf(jj,kk,ll) : 0.0_rt;
178  };
179 
180  int jg = amrex::coarsen(j,2);
181  Real wx = (j == jg*2) ? 0.0_rt : 0.5_rt;
182  Real owx = 1.0_rt-wx;
183 
184  int kg = amrex::coarsen(k,2);
185  Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
186  Real owy = 1.0_rt-wy;
187 
188 #if defined(WARPX_DIM_1D_Z)
189 
190  // interp from coarse nodal to fine nodal
191  Real bg = owx * Byg(jg ,0,0)
192  + wx * Byg(jg+1,0,0);
193 
194  // interp from coarse staggered to fine nodal
195  Real bc = owx * Byc(jg ,0,0)
196  + wx * Byc(jg+1,0,0);
197 
198  // interp from fine staggered to fine nodal
199  Real bf = 0.5_rt*(Byf_zeropad(j-1,0,0) + Byf_zeropad(j,0,0));
200  amrex::ignore_unused(owy);
201 
202 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
203 
204  // interp from coarse nodal to fine nodal
205  Real bg = owx * owy * Byg(jg ,kg ,0)
206  + owx * wy * Byg(jg ,kg+1,0)
207  + wx * owy * Byg(jg+1,kg ,0)
208  + wx * wy * Byg(jg+1,kg+1,0);
209 
210  // interp from coarse stagged (cell-centered for By) to fine nodal
211  wx = 0.5_rt-wx; owx = 1.0_rt-wx;
212  wy = 0.5_rt-wy; owy = 1.0_rt-wy;
213  Real bc = owx * owy * Byc(jg ,kg ,0)
214  + owx * wy * Byc(jg ,kg-1,0)
215  + wx * owy * Byc(jg-1,kg ,0)
216  + wx * wy * Byc(jg-1,kg-1,0);
217 
218  // interp form fine stagger (cell-centered for By) to fine nodal
219  Real bf = 0.25_rt*(Byf_zeropad(j,k,0) + Byf_zeropad(j-1,k,0) + Byf_zeropad(j,k-1,0) + Byf_zeropad(j-1,k-1,0));
220 
221 #else
222 
223  int lg = amrex::coarsen(l,2);
224  Real wz = (l == lg*2) ? 0.0_rt : 0.5_rt;
225  Real owz = 1.0_rt-wz;
226 
227  // interp from coarse nodal to fine nodal
228  Real bg = owx * owy * owz * Byg(jg ,kg ,lg )
229  + wx * owy * owz * Byg(jg+1,kg ,lg )
230  + owx * wy * owz * Byg(jg ,kg+1,lg )
231  + wx * wy * owz * Byg(jg+1,kg+1,lg )
232  + owx * owy * wz * Byg(jg ,kg ,lg+1)
233  + wx * owy * wz * Byg(jg+1,kg ,lg+1)
234  + owx * wy * wz * Byg(jg ,kg+1,lg+1)
235  + wx * wy * wz * Byg(jg+1,kg+1,lg+1);
236 
237  // interp from coarse staggered to fine nodal
238  wx = 0.5_rt-wx; owx = 1.0_rt-wx;
239  wz = 0.5_rt-wz; owz = 1.0_rt-wz;
240  Real bc = owx * owy * owz * Byc(jg ,kg ,lg )
241  + wx * owy * owz * Byc(jg-1,kg ,lg )
242  + owx * wy * owz * Byc(jg ,kg+1,lg )
243  + wx * wy * owz * Byc(jg-1,kg+1,lg )
244  + owx * owy * wz * Byc(jg ,kg ,lg-1)
245  + wx * owy * wz * Byc(jg-1,kg ,lg-1)
246  + owx * wy * wz * Byc(jg ,kg+1,lg-1)
247  + wx * wy * wz * Byc(jg-1,kg+1,lg-1);
248 
249  // interp from fine stagged to fine nodal
250  Real bf = 0.25_rt*(Byf_zeropad(j-1,k,l-1) + Byf_zeropad(j,k,l-1) + Byf_zeropad(j-1,k,l) + Byf_zeropad(j,k,l));
251 
252 #endif
253 
254  Bya(j,k,l) = bg + (bf-bc);
255 }
256 
257 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
258 void warpx_interp_nd_bfield_z (int j, int k, int l,
259  amrex::Array4<amrex::Real> const& Bza,
260  amrex::Array4<amrex::Real const> const& Bzf,
261  amrex::Array4<amrex::Real const> const& Bzc,
262  amrex::Array4<amrex::Real const> const& Bzg)
263 {
264  using namespace amrex;
265 
266  // Pad Bzf with zeros beyond ghost cells for out-of-bound accesses
267  const auto Bzf_zeropad = [Bzf] (const int jj, const int kk, const int ll) noexcept
268  {
269  return Bzf.contains(jj,kk,ll) ? Bzf(jj,kk,ll) : 0.0_rt;
270  };
271 
272  int jg = amrex::coarsen(j,2);
273  Real wx = (j == jg*2) ? 0.0_rt : 0.5_rt;
274  Real owx = 1.0_rt-wx;
275 
276  int kg = amrex::coarsen(k,2);
277  Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
278  Real owy = 1.0_rt-wy;
279 
280 #if defined(WARPX_DIM_1D_Z)
281 
282  // interp from coarse nodal to fine nodal
283  Real bg = owx * Bzg(jg ,0,0)
284  + wx * Bzg(jg+1,0,0);
285 
286  // interp from coarse staggered to fine nodal
287  Real bc = owx * Bzc(jg ,0,0)
288  + wx * Bzc(jg+1,0,0);
289 
290  // interp from fine staggered to fine nodal
291  Real bf = Bzf_zeropad(j,0,0);
292  amrex::ignore_unused(owy);
293 
294 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
295 
296  // interp from coarse nodal to fine nodal
297  Real bg = owx * owy * Bzg(jg ,kg ,0)
298  + owx * wy * Bzg(jg ,kg+1,0)
299  + wx * owy * Bzg(jg+1,kg ,0)
300  + wx * wy * Bzg(jg+1,kg+1,0);
301 
302  // interp from coarse staggered to fine nodal
303  wx = 0.5_rt-wx; owx = 1.0_rt-wx;
304  Real bc = owx * owy * Bzc(jg ,kg ,0)
305  + owx * wy * Bzc(jg ,kg+1,0)
306  + wx * owy * Bzc(jg-1,kg ,0)
307  + wx * wy * Bzc(jg-1,kg+1,0);
308 
309  // interp from fine staggered to fine nodal
310  Real bf = 0.5_rt*(Bzf_zeropad(j-1,k,0) + Bzf_zeropad(j,k,0));
311 
312 #else
313 
314  int lg = amrex::coarsen(l,2);
315  Real wz = (l == lg*2) ? 0.0_rt : 0.5_rt;
316  Real owz = 1.0_rt-wz;
317 
318  // interp from coarse nodal to fine nodal
319  Real bg = owx * owy * owz * Bzg(jg ,kg ,lg )
320  + wx * owy * owz * Bzg(jg+1,kg ,lg )
321  + owx * wy * owz * Bzg(jg ,kg+1,lg )
322  + wx * wy * owz * Bzg(jg+1,kg+1,lg )
323  + owx * owy * wz * Bzg(jg ,kg ,lg+1)
324  + wx * owy * wz * Bzg(jg+1,kg ,lg+1)
325  + owx * wy * wz * Bzg(jg ,kg+1,lg+1)
326  + wx * wy * wz * Bzg(jg+1,kg+1,lg+1);
327 
328  // interp from coarse staggered to fine nodal
329  wx = 0.5_rt-wx; owx = 1.0_rt-wx;
330  wy = 0.5_rt-wy; owy = 1.0_rt-wy;
331  Real bc = owx * owy * owz * Bzc(jg ,kg ,lg )
332  + wx * owy * owz * Bzc(jg-1,kg ,lg )
333  + owx * wy * owz * Bzc(jg ,kg-1,lg )
334  + wx * wy * owz * Bzc(jg-1,kg-1,lg )
335  + owx * owy * wz * Bzc(jg ,kg ,lg+1)
336  + wx * owy * wz * Bzc(jg-1,kg ,lg+1)
337  + owx * wy * wz * Bzc(jg ,kg-1,lg+1)
338  + wx * wy * wz * Bzc(jg-1,kg-1,lg+1);
339 
340  // interp from fine stagged to fine nodal
341  Real bf = 0.25_rt*(Bzf_zeropad(j-1,k-1,l) + Bzf_zeropad(j,k-1,l) + Bzf_zeropad(j-1,k,l) + Bzf_zeropad(j,k,l));
342 
343 #endif
344 
345  Bza(j,k,l) = bg + (bf-bc);
346 }
347 
348 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
349 void warpx_interp_nd_efield_x (int j, int k, int l,
350  amrex::Array4<amrex::Real> const& Exa,
351  amrex::Array4<amrex::Real const> const& Exf,
352  amrex::Array4<amrex::Real const> const& Exc,
353  amrex::Array4<amrex::Real const> const& Exg)
354 {
355  using namespace amrex;
356 
357  // Pad Exf with zeros beyond ghost cells for out-of-bound accesses
358  const auto Exf_zeropad = [Exf] (const int jj, const int kk, const int ll) noexcept
359  {
360  return Exf.contains(jj,kk,ll) ? Exf(jj,kk,ll) : 0.0_rt;
361  };
362 
363  int jg = amrex::coarsen(j,2);
364  Real wx = (j == jg*2) ? 0.0_rt : 0.5_rt;
365  Real owx = 1.0_rt-wx;
366 
367  int kg = amrex::coarsen(k,2);
368  Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
369  Real owy = 1.0_rt-wy;
370 
371 #if defined(WARPX_DIM_1D_Z)
372 
373  // interp from coarse nodal to fine nodal
374  Real eg = owx * Exg(jg ,0,0)
375  + wx * Exg(jg+1,0,0);
376 
377  // interp from coarse staggered to fine nodal
378  Real ec = owx * Exc(jg ,0,0)
379  + wx * Exc(jg+1,0,0);
380 
381  // interp from fine staggered to fine nodal
382  Real ef = Exf_zeropad(j,0,0);
383  amrex::ignore_unused(owy);
384 
385 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
386 
387  // interp from coarse nodal to fine nodal
388  Real eg = owx * owy * Exg(jg ,kg ,0)
389  + owx * wy * Exg(jg ,kg+1,0)
390  + wx * owy * Exg(jg+1,kg ,0)
391  + wx * wy * Exg(jg+1,kg+1,0);
392 
393  // interp from coarse staggered to fine nodal
394  wx = 0.5_rt-wx; owx = 1.0_rt-wx;
395  Real ec = owx * owy * Exc(jg ,kg ,0)
396  + owx * wy * Exc(jg ,kg+1,0)
397  + wx * owy * Exc(jg-1,kg ,0)
398  + wx * wy * Exc(jg-1,kg+1,0);
399 
400  // interp from fine staggered to fine nodal
401  Real ef = 0.5_rt*(Exf_zeropad(j-1,k,0) + Exf_zeropad(j,k,0));
402 
403 #else
404 
405  int lg = amrex::coarsen(l,2);
406  Real wz = (l == lg*2) ? 0.0 : 0.5;
407  Real owz = 1.0_rt-wz;
408 
409  // interp from coarse nodal to fine nodal
410  Real eg = owx * owy * owz * Exg(jg ,kg ,lg )
411  + wx * owy * owz * Exg(jg+1,kg ,lg )
412  + owx * wy * owz * Exg(jg ,kg+1,lg )
413  + wx * wy * owz * Exg(jg+1,kg+1,lg )
414  + owx * owy * wz * Exg(jg ,kg ,lg+1)
415  + wx * owy * wz * Exg(jg+1,kg ,lg+1)
416  + owx * wy * wz * Exg(jg ,kg+1,lg+1)
417  + wx * wy * wz * Exg(jg+1,kg+1,lg+1);
418 
419  // interp from coarse staggered to fine nodal
420  wx = 0.5_rt-wx; owx = 1.0_rt-wx;
421  Real ec = owx * owy * owz * Exc(jg ,kg ,lg )
422  + wx * owy * owz * Exc(jg-1,kg ,lg )
423  + owx * wy * owz * Exc(jg ,kg+1,lg )
424  + wx * wy * owz * Exc(jg-1,kg+1,lg )
425  + owx * owy * wz * Exc(jg ,kg ,lg+1)
426  + wx * owy * wz * Exc(jg-1,kg ,lg+1)
427  + owx * wy * wz * Exc(jg ,kg+1,lg+1)
428  + wx * wy * wz * Exc(jg-1,kg+1,lg+1);
429 
430  // interp from fine staggered to fine nodal
431  Real ef = 0.5_rt*(Exf_zeropad(j-1,k,l) + Exf_zeropad(j,k,l));
432 
433 #endif
434 
435  Exa(j,k,l) = eg + (ef-ec);
436 }
437 
438 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
439 void warpx_interp_nd_efield_y (int j, int k, int l,
440  amrex::Array4<amrex::Real> const& Eya,
441  amrex::Array4<amrex::Real const> const& Eyf,
442  amrex::Array4<amrex::Real const> const& Eyc,
443  amrex::Array4<amrex::Real const> const& Eyg)
444 {
445  using namespace amrex;
446 
447  // Pad Eyf with zeros beyond ghost cells for out-of-bound accesses
448  const auto Eyf_zeropad = [Eyf] (const int jj, const int kk, const int ll) noexcept
449  {
450  return Eyf.contains(jj,kk,ll) ? Eyf(jj,kk,ll) : 0.0_rt;
451  };
452 
453  int jg = amrex::coarsen(j,2);
454  Real wx = (j == jg*2) ? 0.0_rt : 0.5_rt;
455  Real owx = 1.0_rt-wx;
456 
457  int kg = amrex::coarsen(k,2);
458  Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
459  Real owy = 1.0_rt-wy;
460 
461 
462 
463 #if defined(WARPX_DIM_1D_Z)
464 
465  // interp from coarse nodal to fine nodal
466  Real eg = owx * Eyg(jg ,0,0)
467  + wx * Eyg(jg+1,0,0);
468 
469  // interp from coarse staggered to fine nodal
470  Real ec = owx * Eyc(jg ,0,0)
471  + wx * Eyc(jg+1,0,0);
472 
473  // interp from fine staggered to fine nodal
474  Real ef = Eyf_zeropad(j,0,0);
475  amrex::ignore_unused(owy);
476 
477 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
478 
479  // interp from coarse nodal to fine nodal
480  Real eg = owx * owy * Eyg(jg ,kg ,0)
481  + owx * wy * Eyg(jg ,kg+1,0)
482  + wx * owy * Eyg(jg+1,kg ,0)
483  + wx * wy * Eyg(jg+1,kg+1,0);
484 
485  // interp from coarse staggered to fine nodal (Eyc is actually nodal)
486  Real ec = owx * owy * Eyc(jg ,kg ,0)
487  + owx * wy * Eyc(jg ,kg+1,0)
488  + wx * owy * Eyc(jg+1,kg ,0)
489  + wx * wy * Eyc(jg+1,kg+1,0);
490 
491  // interp from fine staggered to fine nodal
492  Real ef = Eyf_zeropad(j,k,0);
493 
494 #else
495 
496  int lg = amrex::coarsen(l,2);
497  Real wz = (l == lg*2) ? 0.0 : 0.5;
498  Real owz = 1.0_rt-wz;
499 
500  // interp from coarse nodal to fine nodal
501  Real eg = owx * owy * owz * Eyg(jg ,kg ,lg )
502  + wx * owy * owz * Eyg(jg+1,kg ,lg )
503  + owx * wy * owz * Eyg(jg ,kg+1,lg )
504  + wx * wy * owz * Eyg(jg+1,kg+1,lg )
505  + owx * owy * wz * Eyg(jg ,kg ,lg+1)
506  + wx * owy * wz * Eyg(jg+1,kg ,lg+1)
507  + owx * wy * wz * Eyg(jg ,kg+1,lg+1)
508  + wx * wy * wz * Eyg(jg+1,kg+1,lg+1);
509 
510  // interp from coarse staggered to fine nodal
511  wy = 0.5_rt-wy; owy = 1.0_rt-wy;
512  Real ec = owx * owy * owz * Eyc(jg ,kg ,lg )
513  + wx * owy * owz * Eyc(jg+1,kg ,lg )
514  + owx * wy * owz * Eyc(jg ,kg-1,lg )
515  + wx * wy * owz * Eyc(jg+1,kg-1,lg )
516  + owx * owy * wz * Eyc(jg ,kg ,lg+1)
517  + wx * owy * wz * Eyc(jg+1,kg ,lg+1)
518  + owx * wy * wz * Eyc(jg ,kg-1,lg+1)
519  + wx * wy * wz * Eyc(jg+1,kg-1,lg+1);
520 
521  // interp from fine staggered to fine nodal
522  Real ef = 0.5_rt*(Eyf_zeropad(j,k-1,l) + Eyf_zeropad(j,k,l));
523 
524 #endif
525 
526  Eya(j,k,l) = eg + (ef-ec);
527 }
528 
529 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
530 void warpx_interp_nd_efield_z (int j, int k, int l,
531  amrex::Array4<amrex::Real> const& Eza,
532  amrex::Array4<amrex::Real const> const& Ezf,
533  amrex::Array4<amrex::Real const> const& Ezc,
534  amrex::Array4<amrex::Real const> const& Ezg)
535 {
536  using namespace amrex;
537 
538  // Pad Ezf with zeros beyond ghost cells for out-of-bound accesses
539  const auto Ezf_zeropad = [Ezf] (const int jj, const int kk, const int ll) noexcept
540  {
541  return Ezf.contains(jj,kk,ll) ? Ezf(jj,kk,ll) : 0.0_rt;
542  };
543 
544  int jg = amrex::coarsen(j,2);
545  Real wx = (j == jg*2) ? 0.0_rt : 0.5_rt;
546  Real owx = 1.0_rt-wx;
547 
548  int kg = amrex::coarsen(k,2);
549  Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt;
550  Real owy = 1.0_rt-wy;
551 
552 #if defined(WARPX_DIM_1D_Z)
553 
554  // interp from coarse nodal to fine nodal
555  Real eg = owx * Ezg(jg ,0,0)
556  + wx * Ezg(jg+1,0,0);
557 
558  // interp from coarse staggered to fine nodal
559  Real ec = owx * Ezc(jg ,0,0)
560  + wx * Ezc(jg+1,0,0);
561 
562  // interp from fine staggered to fine nodal
563  Real ef = 0.5_rt*(Ezf_zeropad(j-1,0,0) + Ezf_zeropad(j,0,0));
564  amrex::ignore_unused(owy);
565 
566 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
567 
568  // interp from coarse nodal to fine nodal
569  Real eg = owx * owy * Ezg(jg ,kg ,0)
570  + owx * wy * Ezg(jg ,kg+1,0)
571  + wx * owy * Ezg(jg+1,kg ,0)
572  + wx * wy * Ezg(jg+1,kg+1,0);
573 
574  // interp from coarse stagged to fine nodal
575  wy = 0.5_rt-wy; owy = 1.0_rt-wy;
576  Real ec = owx * owy * Ezc(jg ,kg ,0)
577  + owx * wy * Ezc(jg ,kg-1,0)
578  + wx * owy * Ezc(jg+1,kg ,0)
579  + wx * wy * Ezc(jg+1,kg-1,0);
580 
581  // interp from fine staggered to fine nodal
582  Real ef = 0.5_rt*(Ezf_zeropad(j,k-1,0) + Ezf_zeropad(j,k,0));
583 
584 #else
585 
586  int lg = amrex::coarsen(l,2);
587  Real wz = (l == lg*2) ? 0.0_rt : 0.5_rt;
588  Real owz = 1.0_rt-wz;
589 
590  // interp from coarse nodal to fine nodal
591  Real eg = owx * owy * owz * Ezg(jg ,kg ,lg )
592  + wx * owy * owz * Ezg(jg+1,kg ,lg )
593  + owx * wy * owz * Ezg(jg ,kg+1,lg )
594  + wx * wy * owz * Ezg(jg+1,kg+1,lg )
595  + owx * owy * wz * Ezg(jg ,kg ,lg+1)
596  + wx * owy * wz * Ezg(jg+1,kg ,lg+1)
597  + owx * wy * wz * Ezg(jg ,kg+1,lg+1)
598  + wx * wy * wz * Ezg(jg+1,kg+1,lg+1);
599 
600  // interp from coarse staggered to fine nodal
601  wz = 0.5_rt-wz; owz = 1.0_rt-wz;
602  Real ec = owx * owy * owz * Ezc(jg ,kg ,lg )
603  + wx * owy * owz * Ezc(jg+1,kg ,lg )
604  + owx * wy * owz * Ezc(jg ,kg+1,lg )
605  + wx * wy * owz * Ezc(jg+1,kg+1,lg )
606  + owx * owy * wz * Ezc(jg ,kg ,lg-1)
607  + wx * owy * wz * Ezc(jg+1,kg ,lg-1)
608  + owx * wy * wz * Ezc(jg ,kg+1,lg-1)
609  + wx * wy * wz * Ezc(jg+1,kg+1,lg-1);
610 
611  // interp from fine staggered to fine nodal
612  Real ef = 0.5_rt*(Ezf_zeropad(j,k,l-1) + Ezf_zeropad(j,k,l));
613 
614 #endif
615 
616  Eza(j,k,l) = eg + (ef-ec);
617 }
618 
638 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
639 void warpx_interp (const int j,
640  const int k,
641  const int l,
642  amrex::Array4<amrex::Real > const& dst_arr,
643  amrex::Array4<amrex::Real const> const& src_arr,
644  const amrex::IntVect& dst_stag,
645  const amrex::IntVect& src_stag,
646  const int nox = 2,
647  const int noy = 2,
648  const int noz = 2,
649  amrex::Real const* stencil_coeffs_x = nullptr,
650  amrex::Real const* stencil_coeffs_y = nullptr,
651  amrex::Real const* stencil_coeffs_z = nullptr)
652 {
653  using namespace amrex;
654 
655  // Pad input array with zeros beyond ghost cells
656  // for out-of-bound accesses due to large-stencil operations
657  const auto src_arr_zeropad = [src_arr] (const int jj, const int kk, const int ll) noexcept
658  {
659  return src_arr.contains(jj,kk,ll) ? src_arr(jj,kk,ll) : 0.0_rt;
660  };
661 
662  // Avoid compiler warnings
663 #if defined(WARPX_DIM_1D_Z)
664  amrex::ignore_unused(nox, noy, stencil_coeffs_x, stencil_coeffs_y);
665 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
666  amrex::ignore_unused(noy, stencil_coeffs_y);
667 #endif
668 
669  // If dst_nodal = true , we are centering from a staggered grid to a nodal grid
670  // If dst_nodal = false, we are centering from a nodal grid to a staggered grid
671  const bool dst_nodal = (dst_stag == amrex::IntVect::TheNodeVector());
672 
673  // See 1D examples below to understand the meaning of this integer shift
674  const int shift = (dst_nodal) ? 0 : 1;
675 
676  // Staggering (s = 0 if cell-centered, s = 1 if nodal)
677  const int sj = (dst_nodal) ? src_stag[0] : dst_stag[0];
678 #if (AMREX_SPACEDIM >= 2)
679  const int sk = (dst_nodal) ? src_stag[1] : dst_stag[1];
680 #endif
681 #if defined(WARPX_DIM_3D)
682  const int sl = (dst_nodal) ? src_stag[2] : dst_stag[2];
683 #endif
684 
685  // Interpolate along j,k,l only if source MultiFab is staggered along j,k,l
686  const bool interp_j = (sj == 0);
687 #if (AMREX_SPACEDIM >= 2)
688  const bool interp_k = (sk == 0);
689 #endif
690 #if defined(WARPX_DIM_3D)
691  const bool interp_l = (sl == 0);
692 #endif
693 
694 #if defined(WARPX_DIM_1D_Z)
695  const int noj = noz;
696 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
697  const int noj = nox;
698  const int nok = noz;
699 #elif defined(WARPX_DIM_3D)
700  const int noj = nox;
701  const int nok = noy;
702  const int nol = noz;
703 #endif
704 
705  // Additional normalization factor
706  const amrex::Real wj = (interp_j) ? 0.5_rt : 1.0_rt;
707 #if defined(WARPX_DIM_1D_Z)
708  constexpr amrex::Real wk = 1.0_rt;
709  constexpr amrex::Real wl = 1.0_rt;
710 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
711  const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt;
712  constexpr amrex::Real wl = 1.0_rt;
713 #elif defined(WARPX_DIM_3D)
714  const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt;
715  const amrex::Real wl = (interp_l) ? 0.5_rt : 1.0_rt;
716 #endif
717 
718  // Min and max for interpolation loop along j
719  const int jmin = (interp_j) ? j - noj/2 + shift : j;
720  const int jmax = (interp_j) ? j + noj/2 + shift - 1 : j;
721 
722  // Min and max for interpolation loop along k
723 #if defined(WARPX_DIM_1D_Z)
724  // k = 0 always
725  const int kmin = k;
726  const int kmax = k;
727 #else
728  const int kmin = (interp_k) ? k - nok/2 + shift : k;
729  const int kmax = (interp_k) ? k + nok/2 + shift - 1 : k;
730 #endif
731 
732  // Min and max for interpolation loop along l
733 #if (AMREX_SPACEDIM <= 2)
734  // l = 0 always
735  const int lmin = l;
736  const int lmax = l;
737 #elif defined(WARPX_DIM_3D)
738  const int lmin = (interp_l) ? l - nol/2 + shift : l;
739  const int lmax = (interp_l) ? l + nol/2 + shift - 1 : l;
740 #endif
741 
742  // Number of interpolation points
743  const int nj = jmax - jmin;
744  const int nk = kmax - kmin;
745  const int nl = lmax - lmin;
746 
747  // Example of 1D centering from nodal grid to nodal grid (simple copy):
748  //
749  // j
750  // --o-----o-----o-- result(j) = f(j)
751  // --o-----o-----o--
752  // j-1 j j+1
753  //
754  // Example of 1D linear centering from staggered grid to nodal grid:
755  //
756  // j
757  // --o-----o-----o-- result(j) = (f(j-1) + f(j)) / 2
758  // -----x-----x-----
759  // j-1 j
760  //
761  // Example of 1D linear centering from nodal grid to staggered grid:
762  // (note the shift of +1 in the indices with respect to the case above, see variable "shift")
763  //
764  // j
765  // --x-----x-----x-- result(j) = (f(j) + f(j+1)) / 2
766  // -----o-----o-----
767  // j j+1
768  //
769  // Example of 1D finite-order centering from staggered grid to nodal grid:
770  //
771  // j
772  // --o-----o-----o-----o-----o-----o-----o-- result(j) = c_0 * (f(j-1) + f(j) ) / 2
773  // -----x-----x-----x-----x-----x-----x----- + c_1 * (f(j-2) + f(j+1)) / 2
774  // j-3 j-2 j-1 j j+1 j+2 + c_2 * (f(j-3) + f(j+2)) / 2
775  // c_2 c_1 c_0 c_0 c_1 c_2 + ...
776  //
777  // Example of 1D finite-order centering from nodal grid to staggered grid:
778  // (note the shift of +1 in the indices with respect to the case above, see variable "shift")
779  //
780  // j
781  // --x-----x-----x-----x-----x-----x-----x-- result(j) = c_0 * (f(j) + f(j+1)) / 2
782  // -----o-----o-----o-----o-----o-----o----- + c_1 * (f(j-1) + f(j+2)) / 2
783  // j-2 j-1 j j+1 j+2 j+3 + c_2 * (f(j-2) + f(j+3)) / 2
784  // c_2 c_1 c_0 c_0 c_1 c_2 + ...
785 
786  amrex::Real res = 0.0_rt;
787 
788  amrex::Real cj = 1.0_rt;
789  amrex::Real ck = 1.0_rt;
790  amrex::Real cl = 1.0_rt;
791 
792 #if defined(WARPX_DIM_1D_Z)
793  amrex::Real const* scj = stencil_coeffs_z;
794 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
795  amrex::Real const* scj = stencil_coeffs_x;
796  amrex::Real const* sck = stencil_coeffs_z;
797 #elif defined(WARPX_DIM_3D)
798  amrex::Real const* scj = stencil_coeffs_x;
799  amrex::Real const* sck = stencil_coeffs_y;
800  amrex::Real const* scl = stencil_coeffs_z;
801 #endif
802 
803  for (int ll = 0; ll <= nl; ll++)
804  {
805 #if defined(WARPX_DIM_3D)
806  if (interp_l) cl = scl[ll];
807 #endif
808  for (int kk = 0; kk <= nk; kk++)
809  {
810 #if (AMREX_SPACEDIM >= 2)
811  if (interp_k) ck = sck[kk];
812 #endif
813  for (int jj = 0; jj <= nj; jj++)
814  {
815  if (interp_j) cj = scj[jj];
816 
817  res += cj * ck * cl * src_arr_zeropad(jmin+jj,kmin+kk,lmin+ll);
818  }
819  }
820  }
821 
822  dst_arr(j,k,l) = wj * wk * wl * res;
823 }
824 
825 #endif
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
AMReX_FArrayBox.H
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
warpx_interp_nd_efield_x
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_nd_efield_x(int j, int k, int l, amrex::Array4< amrex::Real > const &Exa, amrex::Array4< amrex::Real const > const &Exf, amrex::Array4< amrex::Real const > const &Exc, amrex::Array4< amrex::Real const > const &Exg)
Definition: WarpXComm_K.H:349
warpx_interp_nd_bfield_x
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_nd_bfield_x(int j, int k, int l, amrex::Array4< amrex::Real > const &Bxa, amrex::Array4< amrex::Real const > const &Bxf, amrex::Array4< amrex::Real const > const &Bxc, amrex::Array4< amrex::Real const > const &Bxg)
Definition: WarpXComm_K.H:76
warpx_interp_nd_efield_y
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_nd_efield_y(int j, int k, int l, amrex::Array4< amrex::Real > const &Eya, amrex::Array4< amrex::Real const > const &Eyf, amrex::Array4< amrex::Real const > const &Eyc, amrex::Array4< amrex::Real const > const &Eyg)
Definition: WarpXComm_K.H:439
warpx_interp
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp(int j, int k, int l, amrex::Array4< amrex::Real > const &arr_aux, amrex::Array4< amrex::Real const > const &arr_fine, amrex::Array4< amrex::Real const > const &arr_coarse, const amrex::IntVect &arr_stag, const amrex::IntVect &rr)
Definition: WarpXComm_K.H:14
warpx_interp_nd_bfield_y
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_nd_bfield_y(int j, int k, int l, amrex::Array4< amrex::Real > const &Bya, amrex::Array4< amrex::Real const > const &Byf, amrex::Array4< amrex::Real const > const &Byc, amrex::Array4< amrex::Real const > const &Byg)
Definition: WarpXComm_K.H:166
warpx_interp_nd_bfield_z
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_nd_bfield_z(int j, int k, int l, amrex::Array4< amrex::Real > const &Bza, amrex::Array4< amrex::Real const > const &Bzf, amrex::Array4< amrex::Real const > const &Bzc, amrex::Array4< amrex::Real const > const &Bzg)
Definition: WarpXComm_K.H:258
warpx_interp_nd_efield_z
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_nd_efield_z(int j, int k, int l, amrex::Array4< amrex::Real > const &Eza, amrex::Array4< amrex::Real const > const &Ezf, amrex::Array4< amrex::Real const > const &Ezc, amrex::Array4< amrex::Real const > const &Ezg)
Definition: WarpXComm_K.H:530
amrex::IntVectND< AMREX_SPACEDIM >::TheNodeVector
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr IntVectND< dim > TheNodeVector() noexcept
amrex::coarsen
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
amrex::ignore_unused
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
amrex::shift
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > shift(const BoxND< dim > &b, int dir, int nzones) noexcept
check_interp_points_and_weights.jj
jj
Definition: check_interp_points_and_weights.py:160
stencil.nox
int nox
Definition: stencil.py:442
stencil.noy
int noy
Definition: stencil.py:443
stencil.noz
int noz
Definition: stencil.py:444
amrex::Array4::contains
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool contains(int i, int j, int k) const noexcept
Generated by doxygen 1.9.1