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