WarpX
FieldReduction.H
Go to the documentation of this file.
1 /* Copyright 2021 Neil Zaim
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDREDUCTION_H_
9 #define WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDREDUCTION_H_
10 
11 #include "ReducedDiags.H"
12 #include "FieldSolver/Fields.H"
13 #include "WarpX.H"
14 
15 #include <ablastr/coarsen/sample.H>
16 
17 #include <AMReX_Array.H>
18 #include <AMReX_Box.H>
19 #include <AMReX_Config.H>
20 #include <AMReX_FArrayBox.H>
21 #include <AMReX_FabArray.H>
22 #include <AMReX_Geometry.H>
23 #include <AMReX_GpuControl.H>
24 #include <AMReX_GpuQualifiers.H>
25 #include <AMReX_IndexType.H>
26 #include <AMReX_MFIter.H>
27 #include <AMReX_MultiFab.H>
29 #include <AMReX_Parser.H>
30 #include <AMReX_REAL.H>
31 #include <AMReX_RealBox.H>
32 #include <AMReX_Reduce.H>
33 #include <AMReX_Tuple.H>
34 
35 #include <memory>
36 #include <string>
37 #include <type_traits>
38 #include <vector>
39 
46 {
47 public:
48 
53  FieldReduction(const std::string& rd_name);
54 
61  void ComputeDiags(int step) final;
62 
67  void BackwardCompatibility();
68 
69 private:
72  static constexpr int m_nvars = 12;
73  std::unique_ptr<amrex::Parser> m_parser;
74 
75  // Type of reduction (e.g. Maximum, Minimum or Sum)
77 
78 public:
79 
87  template<typename ReduceOp>
89  {
90  using namespace amrex::literals;
91 
92  // get a reference to WarpX instance
93  auto & warpx = WarpX::GetInstance();
94 
95  constexpr int lev = 0; // This reduced diag currently does not work with mesh refinement
96 
97  amrex::Geometry const & geom = warpx.Geom(lev);
98  const amrex::RealBox& real_box = geom.ProbDomain();
99  const auto dx = geom.CellSizeArray();
100 
101  // get MultiFab data
102  const amrex::MultiFab & Ex = warpx.getField(warpx::fields::FieldType::Efield_aux, lev,0);
103  const amrex::MultiFab & Ey = warpx.getField(warpx::fields::FieldType::Efield_aux, lev,1);
104  const amrex::MultiFab & Ez = warpx.getField(warpx::fields::FieldType::Efield_aux, lev,2);
105  const amrex::MultiFab & Bx = warpx.getField(warpx::fields::FieldType::Bfield_aux, lev,0);
106  const amrex::MultiFab & By = warpx.getField(warpx::fields::FieldType::Bfield_aux, lev,1);
107  const amrex::MultiFab & Bz = warpx.getField(warpx::fields::FieldType::Bfield_aux, lev,2);
108  const amrex::MultiFab & jx = warpx.getField(warpx::fields::FieldType::current_fp, lev,0);
109  const amrex::MultiFab & jy = warpx.getField(warpx::fields::FieldType::current_fp, lev,1);
110  const amrex::MultiFab & jz = warpx.getField(warpx::fields::FieldType::current_fp, lev,2);
111 
112 
113  // General preparation of interpolation and reduction operations
114  const amrex::GpuArray<int,3> cellCenteredtype{0,0,0};
115  const amrex::GpuArray<int,3> reduction_coarsening_ratio{1,1,1};
116  constexpr int reduction_comp = 0;
117 
118  amrex::ReduceOps<ReduceOp> reduce_op;
119  amrex::ReduceData<amrex::Real> reduce_data(reduce_op);
120  using ReduceTuple = typename decltype(reduce_data)::Type;
121 
122  // Prepare interpolation of field components to cell center
123  // The arrays below store the index type (staggering) of each MultiFab, with the third
124  // component set to zero in the two-dimensional case.
125  auto Extype = amrex::GpuArray<int,3>{0,0,0};
126  auto Eytype = amrex::GpuArray<int,3>{0,0,0};
127  auto Eztype = amrex::GpuArray<int,3>{0,0,0};
128  auto Bxtype = amrex::GpuArray<int,3>{0,0,0};
129  auto Bytype = amrex::GpuArray<int,3>{0,0,0};
130  auto Bztype = amrex::GpuArray<int,3>{0,0,0};
131  auto jxtype = amrex::GpuArray<int,3>{0,0,0};
132  auto jytype = amrex::GpuArray<int,3>{0,0,0};
133  auto jztype = amrex::GpuArray<int,3>{0,0,0};
134  for (int i = 0; i < AMREX_SPACEDIM; ++i){
135  Extype[i] = Ex.ixType()[i];
136  Eytype[i] = Ey.ixType()[i];
137  Eztype[i] = Ez.ixType()[i];
138  Bxtype[i] = Bx.ixType()[i];
139  Bytype[i] = By.ixType()[i];
140  Bztype[i] = Bz.ixType()[i];
141  jxtype[i] = jx.ixType()[i];
142  jytype[i] = jy.ixType()[i];
143  jztype[i] = jz.ixType()[i];
144 
145  }
146 
147  // get parser
148  auto reduction_function_parser = m_parser->compile<m_nvars>();
149 
150  // MFIter loop to interpolate fields to cell center and perform reduction
151 #ifdef AMREX_USE_OMP
152 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
153 #endif
154  for ( amrex::MFIter mfi(Ex, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
155  {
156  // Make the box cell centered in preparation for the interpolation (and to avoid
157  // including ghost cells in the calculation)
158  const amrex::Box& box = enclosedCells(mfi.nodaltilebox());
159  const auto& arrEx = Ex[mfi].array();
160  const auto& arrEy = Ey[mfi].array();
161  const auto& arrEz = Ez[mfi].array();
162  const auto& arrBx = Bx[mfi].array();
163  const auto& arrBy = By[mfi].array();
164  const auto& arrBz = Bz[mfi].array();
165  const auto& arrjx = jx[mfi].array();
166  const auto& arrjy = jy[mfi].array();
167  const auto& arrjz = jz[mfi].array();
168 
169  reduce_op.eval(box, reduce_data,
170  [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
171  {
172  // 0.5 is here because position are computed on the cell centers
173 
174 #if defined(WARPX_DIM_1D_Z)
175  const amrex::Real x = 0._rt;
176  const amrex::Real y = 0._rt;
177  const amrex::Real z = (k + 0.5_rt)*dx[0] + real_box.lo(0);
178 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
179  const amrex::Real x = (i + 0.5_rt)*dx[0] + real_box.lo(0);
180  const amrex::Real y = 0._rt;
181  const amrex::Real z = (j + 0.5_rt)*dx[1] + real_box.lo(1);
182 #else
183  const amrex::Real x = (i + 0.5_rt)*dx[0] + real_box.lo(0);
184  const amrex::Real y = (j + 0.5_rt)*dx[1] + real_box.lo(1);
185  const amrex::Real z = (k + 0.5_rt)*dx[2] + real_box.lo(2);
186 #endif
187  const amrex::Real Ex_interp = ablastr::coarsen::sample::Interp(arrEx, Extype, cellCenteredtype,
188  reduction_coarsening_ratio, i, j, k, reduction_comp);
189  const amrex::Real Ey_interp = ablastr::coarsen::sample::Interp(arrEy, Eytype, cellCenteredtype,
190  reduction_coarsening_ratio, i, j, k, reduction_comp);
191  const amrex::Real Ez_interp = ablastr::coarsen::sample::Interp(arrEz, Eztype, cellCenteredtype,
192  reduction_coarsening_ratio, i, j, k, reduction_comp);
193  const amrex::Real Bx_interp = ablastr::coarsen::sample::Interp(arrBx, Bxtype, cellCenteredtype,
194  reduction_coarsening_ratio, i, j, k, reduction_comp);
195  const amrex::Real By_interp = ablastr::coarsen::sample::Interp(arrBy, Bytype, cellCenteredtype,
196  reduction_coarsening_ratio, i, j, k, reduction_comp);
197  const amrex::Real Bz_interp = ablastr::coarsen::sample::Interp(arrBz, Bztype, cellCenteredtype,
198  reduction_coarsening_ratio, i, j, k, reduction_comp);
199  const amrex::Real jx_interp = ablastr::coarsen::sample::Interp(arrjx, jxtype, cellCenteredtype,
200  reduction_coarsening_ratio, i, j, k, reduction_comp);
201  const amrex::Real jy_interp = ablastr::coarsen::sample::Interp(arrjy, jytype, cellCenteredtype,
202  reduction_coarsening_ratio, i, j, k, reduction_comp);
203  const amrex::Real jz_interp = ablastr::coarsen::sample::Interp(arrjz, jztype, cellCenteredtype,
204  reduction_coarsening_ratio, i, j, k, reduction_comp);
205 
206  return reduction_function_parser(x, y, z, Ex_interp, Ey_interp, Ez_interp,
207  Bx_interp, By_interp, Bz_interp,
208  jx_interp, jy_interp, jz_interp);
209  });
210  }
211 
212  amrex::Real reduce_value = amrex::get<0>(reduce_data.value());
213 
214  // MPI reduce
216  {
218  }
220  {
222  }
224  {
226  // If reduction operation is a sum, multiply the value by the cell volume so that the
227  // result is the integral of the function over the simulation domain.
228 #if defined(WARPX_DIM_1D_Z)
229  reduce_value *= dx[0];
230 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
231  reduce_value *= dx[0]*dx[1];
232 #else
233  reduce_value *= dx[0]*dx[1]*dx[2];
234 #endif
235  }
236 
237  // Fill output array
238  m_data[0] = reduce_value;
239 
240  // m_data now contains an up-to-date value of the reduced field quantity
241  }
242 
243 };
244 
245 #endif // WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDREDUCTION_H_
#define AMREX_GPU_DEVICE
Definition: FieldReduction.H:46
static constexpr int m_nvars
Definition: FieldReduction.H:72
FieldReduction(const std::string &rd_name)
Definition: FieldReduction.cpp:25
int m_reduction_type
Definition: FieldReduction.H:76
void BackwardCompatibility()
Definition: FieldReduction.cpp:90
void ComputeFieldReduction()
Definition: FieldReduction.H:88
void ComputeDiags(int step) final
Definition: FieldReduction.cpp:103
std::unique_ptr< amrex::Parser > m_parser
Definition: FieldReduction.H:73
Definition: ReducedDiags.H:24
std::vector< amrex::Real > m_data
output data
Definition: ReducedDiags.H:49
static WarpX & GetInstance()
Definition: WarpX.cpp:239
GpuArray< Real, AMREX_SPACEDIM > CellSizeArray() const noexcept
IndexType ixType() const noexcept
Array4< typename FabArray< FArrayBox >::value_type const > array(const MFIter &mfi) const noexcept
const RealBox & ProbDomain() const noexcept
AMREX_GPU_HOST_DEVICE const Real * lo() const &noexcept
std::enable_if_t< IsFabArray< MF >::value > eval(MF const &mf, IntVect const &nghost, D &reduce_data, F &&f)
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, const int i, const int j, const int k, const int comp)
Interpolates the floating point data contained in the source Array4 arr_src, extracted from a fine Mu...
Definition: sample.H:49
void ReduceRealMin(Vector< std::reference_wrapper< Real > > const &)
void ReduceRealSum(Vector< std::reference_wrapper< Real > > const &)
void ReduceRealMax(Vector< std::reference_wrapper< Real > > const &)
bool TilingIfNotGPU() noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box enclosedCells(const Box &b, int dir) noexcept
i
Definition: check_interp_points_and_weights.py:174
tuple dx
lab frame
Definition: stencil.py:429
value
Definition: updateAMReX.py:141
Definition: Fields.H:11