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 
28 void warpx_interp (int j, int k, int l,
29  amrex::Array4<amrex::Real > const& arr_aux,
30  amrex::Array4<amrex::Real const> const& arr_fine,
31  amrex::Array4<amrex::Real const> const& arr_coarse,
32  const amrex::IntVect& arr_stag,
33  const amrex::IntVect& rr)
34 {
35  using namespace amrex;
36 
37  // Pad arr_coarse with zeros beyond ghost cells for out-of-bound accesses
38  const auto arr_coarse_zeropad = [arr_coarse] (const int jj, const int kk, const int ll) noexcept
39  {
40  return arr_coarse.contains(jj,kk,ll) ? arr_coarse(jj,kk,ll) : 0.0_rt;
41  };
42 
43  // NOTE Indices (j,k,l) in the following refer to (z,-,-) in 1D, (x,z,-) in 2D, and (x,y,z) in 3D
44 
45  // Refinement ratio
46  const int rj = rr[0];
47  const int rk = (AMREX_SPACEDIM == 1) ? 1 : rr[1];
48  const int rl = (AMREX_SPACEDIM <= 2) ? 1 : rr[2];
49 
50  // Staggering (0: cell-centered; 1: nodal)
51  const int sj = arr_stag[0];
52  const int sk = (AMREX_SPACEDIM == 1) ? 0 : arr_stag[1];
53  const int sl = (AMREX_SPACEDIM <= 2) ? 0 : arr_stag[2];
54 
55  // Number of points used for interpolation from coarse grid to fine grid
56  const int nj = 2;
57  const int nk = 2;
58  const int nl = 2;
59 
60  const int jc = (sj == 0) ? amrex::coarsen(j - rj/2, rj) : amrex::coarsen(j, rj);
61  const int kc = (sk == 0) ? amrex::coarsen(k - rk/2, rk) : amrex::coarsen(k, rk);
62  const int lc = (sl == 0) ? amrex::coarsen(l - rl/2, rl) : amrex::coarsen(l, rl);
63 
64  // Interpolate from coarse grid to fine grid using 2 points
65  // with weights depending on the distance, for both nodal and cell-centered grids
66  const amrex::Real hj = (sj == 0) ? 0.5_rt : 0._rt;
67  const amrex::Real hk = (sk == 0) ? 0.5_rt : 0._rt;
68  const amrex::Real hl = (sl == 0) ? 0.5_rt : 0._rt;
69 
70  amrex::Real res = 0.0_rt;
71 
72  for (int jj = 0; jj < nj; jj++) {
73  for (int kk = 0; kk < nk; kk++) {
74  for (int ll = 0; ll < nl; ll++) {
75  const amrex::Real wj = (rj - amrex::Math::abs(j + hj - (jc + jj + hj) * rj)) / static_cast<amrex::Real>(rj);
76  const amrex::Real wk = (rk - amrex::Math::abs(k + hk - (kc + kk + hk) * rk)) / static_cast<amrex::Real>(rk);
77  const amrex::Real wl = (rl - amrex::Math::abs(l + hl - (lc + ll + hl) * rl)) / static_cast<amrex::Real>(rl);
78  res += wj * wk * wl * arr_coarse_zeropad(jc+jj,kc+kk,lc+ll);
79  }
80  }
81  }
82  arr_aux(j,k,l) = arr_fine(j,k,l) + res;
83 }
84 
102 void warpx_interp (int j, int k, int l,
103  amrex::Array4<amrex::Real > const& arr_aux,
104  amrex::Array4<amrex::Real const> const& arr_fine,
105  amrex::Array4<amrex::Real const> const& arr_coarse,
106  amrex::Array4<amrex::Real const> const& arr_tmp,
107  const amrex::IntVect& arr_fine_stag,
108  const amrex::IntVect& arr_coarse_stag,
109  const amrex::IntVect& rr)
110 {
111  using namespace amrex;
112 
113  // Pad input arrays with zeros beyond ghost cells
114  // for out-of-bound accesses due to large-stencil operations
115  const auto arr_fine_zeropad = [arr_fine] (const int jj, const int kk, const int ll) noexcept
116  {
117  return arr_fine.contains(jj,kk,ll) ? arr_fine(jj,kk,ll) : 0.0_rt;
118  };
119  const auto arr_coarse_zeropad = [arr_coarse] (const int jj, const int kk, const int ll) noexcept
120  {
121  return arr_coarse.contains(jj,kk,ll) ? arr_coarse(jj,kk,ll) : 0.0_rt;
122  };
123  const auto arr_tmp_zeropad = [arr_tmp] (const int jj, const int kk, const int ll) noexcept
124  {
125  return arr_tmp.contains(jj,kk,ll) ? arr_tmp(jj,kk,ll) : 0.0_rt;
126  };
127 
128  // NOTE Indices (j,k,l) in the following refer to:
129  // - (z,-,-) in 1D
130  // - (x,z,-) in 2D
131  // - (r,z,-) in RZ
132  // - (x,y,z) in 3D
133 
134  // Refinement ratio
135  const int rj = rr[0];
136 #if defined(WARPX_DIM_1D_Z)
137  constexpr int rk = 1;
138  constexpr int rl = 1;
139 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
140  const int rk = rr[1];
141  constexpr int rl = 1;
142 #else
143  const int rk = rr[1];
144  const int rl = rr[2];
145 #endif
146 
147  // Staggering of fine array (0: cell-centered; 1: nodal)
148  const int sj_fp = arr_fine_stag[0];
149 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
150  const int sk_fp = arr_fine_stag[1];
151 #elif defined(WARPX_DIM_3D)
152  const int sk_fp = arr_fine_stag[1];
153  const int sl_fp = arr_fine_stag[2];
154 #endif
155 
156  // Staggering of coarse array (0: cell-centered; 1: nodal)
157  const int sj_cp = arr_coarse_stag[0];
158 #if defined(WARPX_DIM_1D_Z)
159  constexpr int sk_cp = 0;
160  constexpr int sl_cp = 0;
161 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
162  const int sk_cp = arr_coarse_stag[1];
163  constexpr int sl_cp = 0;
164 #else
165  const int sk_cp = arr_coarse_stag[1];
166  const int sl_cp = arr_coarse_stag[2];
167 #endif
168 
169  // Number of points used for interpolation from coarse grid to fine grid
170  int nj;
171  int nk;
172  int nl;
173 
174  int jc = amrex::coarsen(j, rj);
175  int kc = amrex::coarsen(k, rk);
176  int lc = amrex::coarsen(l, rl);
177 
178  amrex::Real tmp = 0.0_rt;
179  amrex::Real fine = 0.0_rt;
180  amrex::Real coarse = 0.0_rt;
181 
182  // 1) Interpolation from coarse nodal to fine nodal
183 
184  nj = 2;
185 #if defined(WARPX_DIM_1D_Z)
186  nk = 1;
187  nl = 1;
188 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
189  nk = 2;
190  nl = 1;
191 #else
192  nk = 2;
193  nl = 2;
194 #endif
195 
196  for (int jj = 0; jj < nj; jj++) {
197  for (int kk = 0; kk < nk; kk++) {
198  for (int ll = 0; ll < nl; ll++) {
199  auto c = arr_tmp_zeropad(jc+jj,kc+kk,lc+ll);
200  c *= (rj - amrex::Math::abs(j - (jc + jj) * rj)) / static_cast<amrex::Real>(rj);
201 #if (AMREX_SPACEDIM >= 2)
202  c *= (rk - amrex::Math::abs(k - (kc + kk) * rk)) / static_cast<amrex::Real>(rk);
203 #endif
204 #if (AMREX_SPACEDIM == 3)
205  c *= (rl - amrex::Math::abs(l - (lc + ll) * rl)) / static_cast<amrex::Real>(rl);
206 #endif
207  tmp += c;
208  }
209  }
210  }
211 
212  // 2) Interpolation from coarse staggered to fine nodal
213 
214  nj = 2;
215 #if defined(WARPX_DIM_1D_Z)
216  nk = 1;
217  nl = 1;
218 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
219  nk = 2;
220  nl = 1;
221 #else
222  nk = 2;
223  nl = 2;
224 #endif
225 
226  const int jn = (sj_cp == 1) ? j : j - rj / 2;
227  const int kn = (sk_cp == 1) ? k : k - rk / 2;
228  const int ln = (sl_cp == 1) ? l : l - rl / 2;
229 
230  jc = amrex::coarsen(jn, rj);
231  kc = amrex::coarsen(kn, rk);
232  lc = amrex::coarsen(ln, rl);
233 
234  for (int jj = 0; jj < nj; jj++) {
235  for (int kk = 0; kk < nk; kk++) {
236  for (int ll = 0; ll < nl; ll++) {
237  auto c = arr_coarse_zeropad(jc+jj,kc+kk,lc+ll);
238  c *= (rj - amrex::Math::abs(jn - (jc + jj) * rj)) / static_cast<amrex::Real>(rj);
239 #if (AMREX_SPACEDIM >= 2)
240  c *= (rk - amrex::Math::abs(kn - (kc + kk) * rk)) / static_cast<amrex::Real>(rk);
241 #endif
242 #if (AMREX_SPACEDIM == 3)
243  c *= (rl - amrex::Math::abs(ln - (lc + ll) * rl)) / static_cast<amrex::Real>(rl);
244 #endif
245  coarse += c;
246  }
247  }
248  }
249 
250  // 3) Interpolation from fine staggered to fine nodal
251 
252  nj = (sj_fp == 0) ? 2 : 1;
253 #if defined(WARPX_DIM_1D_Z)
254  nk = 1;
255  nl = 1;
256 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
257  nk = (sk_fp == 0) ? 2 : 1;
258  nl = 1;
259 #else
260  nk = (sk_fp == 0) ? 2 : 1;
261  nl = (sl_fp == 0) ? 2 : 1;
262 #endif
263 
264  const int jm = (sj_fp == 0) ? j-1 : j;
265 #if defined(WARPX_DIM_1D_Z)
266  const int km = k;
267  const int lm = l;
268 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
269  const int km = (sk_fp == 0) ? k-1 : k;
270  const int lm = l;
271 #else
272  const int km = (sk_fp == 0) ? k-1 : k;
273  const int lm = (sl_fp == 0) ? l-1 : l;
274 #endif
275 
276  for (int jj = 0; jj < nj; jj++) {
277  for (int kk = 0; kk < nk; kk++) {
278  for (int ll = 0; ll < nl; ll++) {
279  fine += arr_fine_zeropad(jm+jj,km+kk,lm+ll);
280  }
281  }
282  }
283  fine = fine/static_cast<amrex::Real>(nj*nk*nl);
284 
285  // Final result
286  arr_aux(j,k,l) = tmp + (fine - coarse);
287 }
302 void warpx_interp (int j, int k, int l,
303  amrex::Array4<amrex::Real > const& arr_aux,
304  amrex::Array4<amrex::Real const> const& arr_fine,
305  const amrex::IntVect& arr_fine_stag)
306 {
307  using namespace amrex;
308 
309  // Pad input arrays with zeros beyond ghost cells
310  // for out-of-bound accesses due to large-stencil operations
311  const auto arr_fine_zeropad = [arr_fine] (const int jj, const int kk, const int ll) noexcept
312  {
313  return arr_fine.contains(jj,kk,ll) ? arr_fine(jj,kk,ll) : 0.0_rt;
314  };
315 
316  // NOTE Indices (j,k,l) in the following refer to:
317  // - (z,-,-) in 1D
318  // - (x,z,-) in 2D
319  // - (r,z,-) in RZ
320  // - (x,y,z) in 3D
321 
322  // Staggering of fine array (0: cell-centered; 1: nodal)
323  const int sj_fp = arr_fine_stag[0];
324 #if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
325  const int sk_fp = arr_fine_stag[1];
326 #elif defined(WARPX_DIM_3D)
327  const int sk_fp = arr_fine_stag[1];
328  const int sl_fp = arr_fine_stag[2];
329 #endif
330 
331  // Number of points used for interpolation from coarse grid to fine grid
332  int nj;
333  int nk;
334  int nl;
335 
336  amrex::Real fine = 0.0_rt;
337 
338  // 3) Interpolation from fine staggered to fine nodal
339 
340  nj = (sj_fp == 0) ? 2 : 1;
341 #if defined(WARPX_DIM_1D_Z)
342  nk = 1;
343  nl = 1;
344 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
345  nk = (sk_fp == 0) ? 2 : 1;
346  nl = 1;
347 #else
348  nk = (sk_fp == 0) ? 2 : 1;
349  nl = (sl_fp == 0) ? 2 : 1;
350 #endif
351 
352  const int jm = (sj_fp == 0) ? j-1 : j;
353 #if defined(WARPX_DIM_1D_Z)
354  const int km = k;
355  const int lm = l;
356 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
357  const int km = (sk_fp == 0) ? k-1 : k;
358  const int lm = l;
359 #else
360  const int km = (sk_fp == 0) ? k-1 : k;
361  const int lm = (sl_fp == 0) ? l-1 : l;
362 #endif
363 
364  for (int jj = 0; jj < nj; jj++) {
365  for (int kk = 0; kk < nk; kk++) {
366  for (int ll = 0; ll < nl; ll++) {
367  fine += arr_fine_zeropad(jm+jj,km+kk,lm+ll);
368  }
369  }
370  }
371  fine = fine/static_cast<amrex::Real>(nj*nk*nl);
372 
373  // Final result
374  arr_aux(j,k,l) = fine;
375 }
376 
397 void warpx_interp (const int j,
398  const int k,
399  const int l,
400  amrex::Array4<amrex::Real > const& dst_arr,
401  amrex::Array4<amrex::Real const> const& src_arr,
402  const amrex::IntVect& dst_stag,
403  const amrex::IntVect& src_stag,
404  const int nox = 2,
405  const int noy = 2,
406  const int noz = 2,
407  amrex::Real const* stencil_coeffs_x = nullptr,
408  amrex::Real const* stencil_coeffs_y = nullptr,
409  amrex::Real const* stencil_coeffs_z = nullptr)
410 {
411  using namespace amrex;
412 
413  // Pad input array with zeros beyond ghost cells
414  // for out-of-bound accesses due to large-stencil operations
415  const auto src_arr_zeropad = [src_arr] (const int jj, const int kk, const int ll) noexcept
416  {
417  return src_arr.contains(jj,kk,ll) ? src_arr(jj,kk,ll) : 0.0_rt;
418  };
419 
420  // Avoid compiler warnings
421 #if defined(WARPX_DIM_1D_Z)
422  amrex::ignore_unused(nox, noy, stencil_coeffs_x, stencil_coeffs_y);
423 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
424  amrex::ignore_unused(noy, stencil_coeffs_y);
425 #endif
426 
427  // If dst_nodal = true , we are centering from a staggered grid to a nodal grid
428  // If dst_nodal = false, we are centering from a nodal grid to a staggered grid
429  const bool dst_nodal = (dst_stag == amrex::IntVect::TheNodeVector());
430 
431  // See 1D examples below to understand the meaning of this integer shift
432  const int shift = (dst_nodal) ? 0 : 1;
433 
434  // Staggering (s = 0 if cell-centered, s = 1 if nodal)
435  const int sj = (dst_nodal) ? src_stag[0] : dst_stag[0];
436 #if (AMREX_SPACEDIM >= 2)
437  const int sk = (dst_nodal) ? src_stag[1] : dst_stag[1];
438 #endif
439 #if defined(WARPX_DIM_3D)
440  const int sl = (dst_nodal) ? src_stag[2] : dst_stag[2];
441 #endif
442 
443  // Interpolate along j,k,l only if source MultiFab is staggered along j,k,l
444  const bool interp_j = (sj == 0);
445 #if (AMREX_SPACEDIM >= 2)
446  const bool interp_k = (sk == 0);
447 #endif
448 #if defined(WARPX_DIM_3D)
449  const bool interp_l = (sl == 0);
450 #endif
451 
452 #if defined(WARPX_DIM_1D_Z)
453  const int noj = noz;
454 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
455  const int noj = nox;
456  const int nok = noz;
457 #elif defined(WARPX_DIM_3D)
458  const int noj = nox;
459  const int nok = noy;
460  const int nol = noz;
461 #endif
462 
463  // Additional normalization factor
464  const amrex::Real wj = (interp_j) ? 0.5_rt : 1.0_rt;
465 #if defined(WARPX_DIM_1D_Z)
466  constexpr amrex::Real wk = 1.0_rt;
467  constexpr amrex::Real wl = 1.0_rt;
468 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
469  const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt;
470  constexpr amrex::Real wl = 1.0_rt;
471 #elif defined(WARPX_DIM_3D)
472  const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt;
473  const amrex::Real wl = (interp_l) ? 0.5_rt : 1.0_rt;
474 #endif
475 
476  // Min and max for interpolation loop along j
477  const int jmin = (interp_j) ? j - noj/2 + shift : j;
478  const int jmax = (interp_j) ? j + noj/2 + shift - 1 : j;
479 
480  // Min and max for interpolation loop along k
481 #if defined(WARPX_DIM_1D_Z)
482  // k = 0 always
483  const int kmin = k;
484  const int kmax = k;
485 #else
486  const int kmin = (interp_k) ? k - nok/2 + shift : k;
487  const int kmax = (interp_k) ? k + nok/2 + shift - 1 : k;
488 #endif
489 
490  // Min and max for interpolation loop along l
491 #if (AMREX_SPACEDIM <= 2)
492  // l = 0 always
493  const int lmin = l;
494  const int lmax = l;
495 #elif defined(WARPX_DIM_3D)
496  const int lmin = (interp_l) ? l - nol/2 + shift : l;
497  const int lmax = (interp_l) ? l + nol/2 + shift - 1 : l;
498 #endif
499 
500  // Number of interpolation points
501  const int nj = jmax - jmin;
502  const int nk = kmax - kmin;
503  const int nl = lmax - lmin;
504 
505  // Example of 1D centering from nodal grid to nodal grid (simple copy):
506  //
507  // j
508  // --o-----o-----o-- result(j) = f(j)
509  // --o-----o-----o--
510  // j-1 j j+1
511  //
512  // Example of 1D linear centering from staggered grid to nodal grid:
513  //
514  // j
515  // --o-----o-----o-- result(j) = (f(j-1) + f(j)) / 2
516  // -----x-----x-----
517  // j-1 j
518  //
519  // Example of 1D linear centering from nodal grid to staggered grid:
520  // (note the shift of +1 in the indices with respect to the case above, see variable "shift")
521  //
522  // j
523  // --x-----x-----x-- result(j) = (f(j) + f(j+1)) / 2
524  // -----o-----o-----
525  // j j+1
526  //
527  // Example of 1D finite-order centering from staggered grid to nodal grid:
528  //
529  // j
530  // --o-----o-----o-----o-----o-----o-----o-- result(j) = c_0 * (f(j-1) + f(j) ) / 2
531  // -----x-----x-----x-----x-----x-----x----- + c_1 * (f(j-2) + f(j+1)) / 2
532  // j-3 j-2 j-1 j j+1 j+2 + c_2 * (f(j-3) + f(j+2)) / 2
533  // c_2 c_1 c_0 c_0 c_1 c_2 + ...
534  //
535  // Example of 1D finite-order centering from nodal grid to staggered grid:
536  // (note the shift of +1 in the indices with respect to the case above, see variable "shift")
537  //
538  // j
539  // --x-----x-----x-----x-----x-----x-----x-- result(j) = c_0 * (f(j) + f(j+1)) / 2
540  // -----o-----o-----o-----o-----o-----o----- + c_1 * (f(j-1) + f(j+2)) / 2
541  // j-2 j-1 j j+1 j+2 j+3 + c_2 * (f(j-2) + f(j+3)) / 2
542  // c_2 c_1 c_0 c_0 c_1 c_2 + ...
543 
544  amrex::Real res = 0.0_rt;
545 
546 #if defined(WARPX_DIM_1D_Z)
547  amrex::Real const* scj = stencil_coeffs_z;
548 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
549  amrex::Real const* scj = stencil_coeffs_x;
550  amrex::Real const* sck = stencil_coeffs_z;
551 #elif defined(WARPX_DIM_3D)
552  amrex::Real const* scj = stencil_coeffs_x;
553  amrex::Real const* sck = stencil_coeffs_y;
554  amrex::Real const* scl = stencil_coeffs_z;
555 #endif
556 
557  for (int ll = 0; ll <= nl; ll++)
558  {
559 #if defined(WARPX_DIM_3D)
560  const amrex::Real cl = (interp_l)? scl[ll] : 1.0_rt;
561 #else
562  const amrex::Real cl = 1.0_rt;
563 #endif
564  for (int kk = 0; kk <= nk; kk++)
565  {
566 #if (AMREX_SPACEDIM >= 2)
567  const amrex::Real ck = (interp_k)? sck[kk] : 1.0_rt;
568 #else
569  const amrex::Real ck = 1.0_rt;
570 #endif
571  for (int jj = 0; jj <= nj; jj++)
572  {
573  const amrex::Real cj = (interp_j)? scj[jj] : 1.0_rt;
574 
575  res += cj * ck * cl * src_arr_zeropad(jmin+jj,kmin+kk,lmin+ll);
576  }
577  }
578  }
579 
580  dst_arr(j,k,l) = wj * wk * wl * res;
581 }
582 
583 #endif
#define AMREX_FORCE_INLINE
#define AMREX_GPU_DEVICE
Array4< Real > fine
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)
Interpolation function called within WarpX::UpdateAuxilaryDataSameType with electromagnetic solver to...
Definition: WarpXComm_K.H:28
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE IntVectND< dim > TheNodeVector() noexcept
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > shift(const BoxND< dim > &b, int dir, int nzones) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
jj
Definition: check_interp_points_and_weights.py:160
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool contains(int i, int j, int k) const noexcept