WarpX
average.H
Go to the documentation of this file.
1 /* Copyright 2022 Edoardo Zoni, Remi Lehe, Prabhat Kumar, Axel Huebl
2  *
3  * This file is part of ABLASTR.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef ABLASTR_COARSEN_AVERAGE_H_
8 #define ABLASTR_COARSEN_AVERAGE_H_
9 
10 #include <AMReX_Array.H>
11 #include <AMReX_Array4.H>
12 #include <AMReX_BLassert.H>
13 #include <AMReX_Extension.H>
14 #include <AMReX_GpuQualifiers.H>
15 #include <AMReX_Math.H>
16 #include <AMReX_REAL.H>
17 #include <AMReX_BaseFwd.H>
18 
19 #include <cstdlib>
20 
21 
27 {
49  amrex::Real
51  amrex::Array4<amrex::Real const> const &arr_src,
52  amrex::GpuArray<int, 3> const &sf,
53  amrex::GpuArray<int, 3> const &sc,
55  int const i,
56  int const j,
57  int const k,
58  int const comp
59  )
60  {
61  using namespace amrex::literals;
62 
63  AMREX_ASSERT_WITH_MESSAGE(sf[0] == sc[0], "Interp: Staggering for component 0 does not match!");
64  AMREX_ASSERT_WITH_MESSAGE(sf[1] == sc[1], "Interp: Staggering for component 1 does not match!");
65  AMREX_ASSERT_WITH_MESSAGE(sf[2] == sc[2], "Interp: Staggering for component 2 does not match!");
66 
67  // Indices of destination array (coarse)
68  int const ic[3] = {i, j, k};
69 
70  // Number of points and starting indices of source array (fine)
71  int np[3], idx_min[3];
72 
73  // Compute number of points
74  for (int l = 0; l < 3; ++l) {
75  if (cr[l] == 1) {
76  np[l] = 1; // no coarsening
77  } else {
78  np[l] = cr[l] * (1 - sf[l]) * (1 - sc[l]) // cell-centered
79  + (2 * (cr[l] - 1) + 1) * sf[l] * sc[l]; // nodal
80  }
81  }
82 
83  // Compute starting indices of source array (fine)
84  for (int l = 0; l < 3; ++l) {
85  if (cr[l] == 1) {
86  idx_min[l] = ic[l]; // no coarsening
87  } else {
88  idx_min[l] = ic[l] * cr[l] * (1 - sf[l]) * (1 - sc[l]) // cell-centered
89  + (ic[l] * cr[l] - cr[l] + 1) * sf[l] * sc[l]; // nodal
90  }
91  }
92 
93  // Auxiliary integer variables
94  int const numx = np[0];
95  int const numy = np[1];
96  int const numz = np[2];
97  int const imin = idx_min[0];
98  int const jmin = idx_min[1];
99  int const kmin = idx_min[2];
100  int const sfx = sf[0];
101  int const sfy = sf[1];
102  int const sfz = sf[2];
103  int const scx = sc[0];
104  int const scy = sc[1];
105  int const scz = sc[2];
106  int const crx = cr[0];
107  int const cry = cr[1];
108  int const crz = cr[2];
109 
110  // Add neutral elements (=0) beyond guard cells in source array (fine)
111  auto const arr_src_safe = [arr_src]
112  AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept {
113  return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt;
114  };
115 
116  // Interpolate over points computed above. Weights are computed in order
117  // to guarantee total charge conservation for both cell-centered data
118  // (equal weights) and nodal data (weights depend on distance between
119  // points on fine and coarse grids). Terms multiplied by (1-sf)*(1-sc)
120  // are ON for cell-centered data and OFF for nodal data, while terms
121  // multiplied by sf*sc are ON for nodal data and OFF for cell-centered data.
122  // Python script Source/Utils/check_interp_points_and_weights.py can be
123  // used to check interpolation points and weights in 1D.
124  amrex::Real c = 0.0_rt;
125  for (int kref = 0; kref < numz; ++kref) {
126  for (int jref = 0; jref < numy; ++jref) {
127  for (int iref = 0; iref < numx; ++iref) {
128  const int ii = imin + iref;
129  const int jj = jmin + jref;
130  const int kk = kmin + kref;
131  const amrex::Real wx = (1.0_rt / static_cast<amrex::Real>(numx)) * (1 - sfx) * (1 - scx) // if cell-centered
132  + ((amrex::Math::abs(crx - amrex::Math::abs(ii - i * crx))) /
133  static_cast<amrex::Real>(crx * crx)) * sfx * scx; // if nodal
134  const amrex::Real wy = (1.0_rt / static_cast<amrex::Real>(numy)) * (1 - sfy) * (1 - scy) // if cell-centered
135  + ((amrex::Math::abs(cry - amrex::Math::abs(jj - j * cry))) /
136  static_cast<amrex::Real>(cry * cry)) * sfy * scy; // if nodal
137  const amrex::Real wz = (1.0_rt / static_cast<amrex::Real>(numz)) * (1 - sfz) * (1 - scz) // if cell-centered
138  + ((amrex::Math::abs(crz - amrex::Math::abs(kk - k * crz))) /
139  static_cast<amrex::Real>(crz * crz)) * sfz * scz; // if nodal
140  c += wx * wy * wz * arr_src_safe(ii, jj, kk, comp);
141  }
142  }
143  }
144  return c;
145  }
146 
160  void
161  Loop (
162  amrex::MultiFab & mf_dst,
163  amrex::MultiFab const & mf_src,
164  int ncomp,
165  amrex::IntVect ngrow,
166  amrex::IntVect crse_ratio
167  );
168 
179  void
180  Coarsen (
181  amrex::MultiFab & mf_dst,
182  amrex::MultiFab const & mf_src,
183  amrex::IntVect crse_ratio
184  );
185 
186 } // namespace ablastr::coarsen::average
187 
188 #endif // ABLASTR_COARSEN_AVERAGE_H_
#define AMREX_ASSERT_WITH_MESSAGE(EX, MSG)
#define AMREX_FORCE_INLINE
#define AMREX_GPU_DEVICE
Definition: average.cpp:24
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Interp(amrex::Array4< amrex::Real const > const &arr_src, amrex::GpuArray< int, 3 > const &sf, amrex::GpuArray< int, 3 > const &sc, amrex::GpuArray< int, 3 > const &cr, int const i, int const j, int const k, int const comp)
Interpolates the floating point data contained in the source Array4 arr_src, extracted from a fine Mu...
Definition: average.H:50
void Coarsen(amrex::MultiFab &mf_dst, amrex::MultiFab const &mf_src, amrex::IntVect const crse_ratio)
Stores in the coarsened MultiFab mf_dst the values obtained by interpolating the data contained in th...
Definition: average.cpp:67
void Loop(amrex::MultiFab &mf_dst, amrex::MultiFab const &mf_src, int const ncomp, amrex::IntVect const ngrow, amrex::IntVect const crse_ratio)
Loops over the boxes of the coarsened MultiFab mf_dst and fills them by interpolating the data contai...
Definition: average.cpp:26
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
constexpr int iz
constexpr int iy
constexpr int ix
i
Definition: check_interp_points_and_weights.py:174
ii
Definition: check_interp_points_and_weights.py:148
jj
Definition: check_interp_points_and_weights.py:160
imin
Definition: check_interp_points_and_weights.py:124
cr
TODO Coarsening for IO: interpolation points and weights def coarsening_points_and_weights_for_IO( i,...
Definition: check_interp_points_and_weights.py:101
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool contains(int i, int j, int k) const noexcept