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