Hipace
Fields.H
Go to the documentation of this file.
1 /* Copyright 2020-2022
2  *
3  * This file is part of HiPACE++.
4  *
5  * Authors: AlexanderSinn, Axel Huebl, MaxThevenet, Remi Lehe
6  * Severin Diederichs, WeiqunZhang, coulibaly-mouhamed
7  * License: BSD-3-Clause-LBNL
8  */
9 #ifndef FIELDS_H_
10 #define FIELDS_H_
11 
13 #include "diagnostics/Diagnostic.H"
14 #include "laser/MultiLaser.H"
15 #include "utils/GPUUtil.H"
16 
17 #include <AMReX_MultiFab.H>
18 #include <AMReX_Vector.H>
19 #include <AMReX_AmrCore.H>
20 
21 #include <algorithm>
22 
23 class Hipace;
24 class MultiLaser;
25 
27 struct WhichSlice {
29 };
30 
31 struct assert_map : std::map<std::string, int> {
32  int operator[] (const std::string& str) const {
33  if (count(str)==0) {
34  amrex::ErrorStream() << "Component '"+str+"' is not allocated\nFields allocated:";
35  for (auto& [alloc_str, num] : *this) {
36  amrex::ErrorStream() << " '" << alloc_str << "' (" << num << "),";
37  }
38  amrex::ErrorStream() << "\n";
39  AMREX_ALWAYS_ASSERT(count(str)!=0);
40  }
41  return at(str);
42  }
43 
44  template<class...Args>
45  void multi_emplace(int& n, Args...comps) {
46  (emplace(comps, n++),...);
47  }
48 };
49 
52 inline std::array<assert_map, WhichSlice::N> Comps{};
53 
56 inline int N_Comps {0};
57 
59 struct Direction{
60  enum dir{x=0, y, z};
61 };
62 
71 inline amrex::Real
72 GetPosOffset (const int dir, const amrex::Geometry& geom, const amrex::Box& box) {
73  using namespace amrex::literals;
74  // match boxes at center point
75  return 0.5_rt * (geom.ProbLo(dir) + geom.ProbHi(dir)
76  - geom.CellSize(dir) * (box.smallEnd(dir) + box.bigEnd(dir)));
77 }
78 
84 class Fields
85 {
86 public:
88  explicit Fields (const int nlev);
89 
98  void AllocData (
99  int lev, amrex::Geometry const& geom, const amrex::BoxArray& slice_ba,
100  const amrex::DistributionMapping& slice_dm, int bin_size);
101 
109  amrex::MultiFab& getSlices (int lev) {return m_slices[lev]; }
113  const amrex::MultiFab& getSlices (int lev) const {return m_slices[lev]; }
119  amrex::MultiFab getField (const int lev, const int islice, const std::string comp) {
120  return amrex::MultiFab(getSlices(lev), amrex::make_alias, Comps[islice][comp], 1);
121  }
125  amrex::MultiFab getStagingArea (const int lev) {
126  return amrex::MultiFab(m_poisson_solver[lev]->StagingArea(), amrex::make_alias, 0, 1);
127  }
131  void checkInit() {
132  for (auto& slices_lev : m_slices) {
134  slices_lev.nComp() == 0 || slices_lev.nGrowVect() == m_slices_nguards,
135  "m_slices[lev].nGrowVect() must be equal to m_slices_nguards"
136  "for initialized slices");
137  }
138  }
147  void Copy (const int lev, const int i_slice, FieldDiagnosticData& fd,
148  const amrex::Geometry& calc_geom, MultiLaser& multi_laser);
149 
156  void InitializeSlices (int lev, int islice, const amrex::Vector<amrex::Geometry>& geom);
157 
166  void ShiftSlices (int lev);
167 
169  void AddRhoIons (const int lev);
170 
181  void SetBoundaryCondition (amrex::Vector<amrex::Geometry> const& geom, const int lev,
182  const int which_slice, std::string component,
183  amrex::MultiFab&& staging_area);
184 
195  void LevelUpBoundary (amrex::Vector<amrex::Geometry> const& geom, const int lev,
196  const int which_slice, const std::string& component,
197  const amrex::IntVect outer_edge, const amrex::IntVect inner_edge);
198 
206  void LevelUp (amrex::Vector<amrex::Geometry> const& geom, const int lev,
207  const int which_slice, const std::string& component);
208 
217  const int current_N_level);
225  void SolvePoissonEz (amrex::Vector<amrex::Geometry> const& geom, const int current_N_level,
226  const int which_slice = WhichSlice::This);
234  void SolvePoissonBxBy (amrex::Vector<amrex::Geometry> const& geom, const int current_N_level,
235  const int which_slice);
244  void SymmetrizeFields (int field_comp, const int lev, const int symm_x, const int symm_y);
253  void InitialBfieldGuess (const amrex::Real relative_Bfield_error,
254  const amrex::Real predcorr_B_error_tolerance, const int lev);
264  void MixAndShiftBfields (const amrex::Real relative_Bfield_error,
265  const amrex::Real relative_Bfield_error_prev_iter,
266  const amrex::Real predcorr_B_mixing_factor, const int lev);
267 
276  amrex::Real ComputeRelBFieldError (const int which_slice, const int which_slice_iter,
277  const amrex::Vector<amrex::Geometry>& geom, const int current_N_level);
278 
279 
288  void InSituComputeDiags (int step, amrex::Real time, int islice, const amrex::Geometry& geom3D,
289  int max_step, amrex::Real max_time);
290 
298  void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom3D,
299  int max_step, amrex::Real max_time);
300 
308  template<class...Args>
309  void setVal (const amrex::Real val, const int lev, const int islice, Args...comps) {
310  static constexpr int ncomps = sizeof...(comps);
311  const amrex::GpuArray<int, ncomps> c_idx = {Comps[islice][comps]...};
312  amrex::MultiFab& mfab = getSlices(lev);
313 
314 #ifdef AMREX_USE_OMP
315 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
316 #endif
317  for (amrex::MFIter mfi(mfab, DfltMfiTlng); mfi.isValid(); ++mfi) {
318  const Array3<amrex::Real> array = mfab.array(mfi);
319  // only one Kernel for all fields
320  amrex::ParallelFor(mfi.growntilebox(),
321  [=] AMREX_GPU_DEVICE(int i, int j, int)
322  {
323  for (int n=0; n<ncomps; ++n) {
324  array(i,j,c_idx[n]) = val;
325  }
326  });
327  }
328  }
329 
337  template<class...Args>
338  void mult (const amrex::Real val, const int lev, const int islice, Args...comps) {
339  static constexpr int ncomps = sizeof...(comps);
340  const amrex::GpuArray<int, ncomps> c_idx = {Comps[islice][comps]...};
341  amrex::MultiFab& mfab = getSlices(lev);
342 
343 #ifdef AMREX_USE_OMP
344 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
345 #endif
346  for (amrex::MFIter mfi(mfab, DfltMfiTlng); mfi.isValid(); ++mfi) {
347  const Array3<amrex::Real> array = mfab.array(mfi);
348  // only one Kernel for all fields
349  amrex::ParallelFor(mfi.growntilebox(),
350  [=] AMREX_GPU_DEVICE(int i, int j, int)
351  {
352  for (int n=0; n<ncomps; ++n) {
353  array(i,j,c_idx[n]) *= val;
354  }
355  });
356  }
357  }
358 
366  template<class...Args>
367  void shift (const int lev, const int islice_dst, const int islice_src, Args...comps) {
368  static constexpr int ncomps = sizeof...(comps);
369  const amrex::GpuArray<int, ncomps> c_idx_src = {Comps[islice_src][comps]...};
370  const amrex::GpuArray<int, ncomps> c_idx_dst = {Comps[islice_dst][comps]...};
371  amrex::MultiFab& mfab = getSlices(lev);
372 
373 #ifdef AMREX_USE_OMP
374 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
375 #endif
376  for (amrex::MFIter mfi(mfab, DfltMfiTlng); mfi.isValid(); ++mfi) {
377  const Array3<amrex::Real> array = mfab.array(mfi);
378  // only one Kernel for all fields
379  amrex::ParallelFor(mfi.growntilebox(),
380  [=] AMREX_GPU_DEVICE(int i, int j, int)
381  {
382  for (int n=0; n<ncomps; ++n) {
383  array(i,j,c_idx_dst[n]) = array(i,j,c_idx_src[n]);
384  }
385  });
386  }
387  }
388 
398  template<int ncomps>
399  void duplicate (const int lev,
400  const int islice_dst, const char * const (&comps_dst)[ncomps],
401  const int islice_src, const char * const (&comps_src)[ncomps]) {
402  amrex::GpuArray<int, ncomps> c_idx_src = {};
403  amrex::GpuArray<int, ncomps> c_idx_dst = {};
404  for (int i=0; i<ncomps; ++i) {
405  c_idx_src[i] = Comps[islice_src][comps_src[i]];
406  c_idx_dst[i] = Comps[islice_dst][comps_dst[i]];
407  }
408  amrex::MultiFab& mfab = getSlices(lev);
409 
410 #ifdef AMREX_USE_OMP
411 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
412 #endif
413  for (amrex::MFIter mfi(mfab, DfltMfiTlng); mfi.isValid(); ++mfi) {
414  const Array3<amrex::Real> array = mfab.array(mfi);
415  // only one Kernel for all fields
416  amrex::ParallelFor(mfi.growntilebox(),
417  [=] AMREX_GPU_DEVICE(int i, int j, int)
418  {
419  for (int n=0; n<ncomps; ++n) {
420  array(i,j,c_idx_dst[n]) = array(i,j,c_idx_src[n]);
421  }
422  });
423  }
424  }
425 
435  template<int ncomps>
436  void add (const int lev,
437  const int islice_dst, const char * const (&comps_dst)[ncomps],
438  const int islice_src, const char * const (&comps_src)[ncomps]) {
439  amrex::GpuArray<int, ncomps> c_idx_src = {};
440  amrex::GpuArray<int, ncomps> c_idx_dst = {};
441  for (int i=0; i<ncomps; ++i) {
442  c_idx_src[i] = Comps[islice_src][comps_src[i]];
443  c_idx_dst[i] = Comps[islice_dst][comps_dst[i]];
444  }
445  amrex::MultiFab& mfab = getSlices(lev);
446 
447 #ifdef AMREX_USE_OMP
448 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
449 #endif
450  for (amrex::MFIter mfi(mfab, DfltMfiTlng); mfi.isValid(); ++mfi) {
451  const Array3<amrex::Real> array = mfab.array(mfi);
452  // only one Kernel for all fields
453  amrex::ParallelFor(mfi.growntilebox(),
454  [=] AMREX_GPU_DEVICE(int i, int j, int)
455  {
456  for (int n=0; n<ncomps; ++n) {
457  array(i,j,c_idx_dst[n]) += array(i,j,c_idx_src[n]);
458  }
459  });
460  }
461  }
462 
470  template<class...Args>
471  void FillBoundary (const amrex::Periodicity& period, const int lev, const int islice,
472  Args...comps) {
473  std::array<int, sizeof...(comps)> c_idx = {Comps[islice][comps]...};
474  amrex::MultiFab& mfab = getSlices(lev);
475 
476  // optimize adjacent fields to one FillBoundary call
477  std::sort(c_idx.begin(), c_idx.end());
478  int scomp = 0;
479  int ncomp = 0;
480  for (unsigned int i=0; i < c_idx.size(); ++i) {
481  if (ncomp==0) {
482  scomp = c_idx[i];
483  ncomp = 1;
484  }
485  if (i+1 >= c_idx.size() || c_idx[i+1] > scomp+ncomp) {
486  mfab.FillBoundary(scomp, ncomp, period);
487  ncomp = 0;
488  } else if (c_idx[i+1] == scomp+ncomp) {
489  ++ncomp;
490  }
491  }
492  }
493 
495  inline static amrex::IntVect m_slices_nguards {-1, -1, -1};
497  inline static amrex::IntVect m_poisson_nguards {-1, -1, -1};
499  bool m_extended_solve = false;
501  bool m_do_symmetrize = false;
502 
503 private:
519  bool m_open_boundary = false;
521  bool m_explicit = false;
525  static constexpr int m_insitu_nrp = 10;
534  std::string m_insitu_file_prefix = "diags/field_insitu";
535 };
536 
537 #endif
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
#define AMREX_ALWAYS_ASSERT(EX)
#define AMREX_GPU_DEVICE
std::array< assert_map, WhichSlice::N > Comps
Map names and indices of each fields in each slice.
Definition: Fields.H:52
amrex::Real GetPosOffset(const int dir, const amrex::Geometry &geom, const amrex::Box &box)
Function to help converting field indexes to positions and backwards. Usage: x = i * dx + GetPosOffse...
Definition: Fields.H:72
int N_Comps
number of fields in each slice
Definition: Fields.H:56
amrex::MFItInfo DfltMfiTlng
Definition: GPUUtil.H:184
Main class handling all field data structures and operations.
Definition: Fields.H:85
void SolvePoissonBxBy(amrex::Vector< amrex::Geometry > const &geom, const int current_N_level, const int which_slice)
Compute Bx and By on the slice container from J by solving two Poisson equations. This function does ...
Definition: Fields.cpp:1012
amrex::Real ComputeRelBFieldError(const int which_slice, const int which_slice_iter, const amrex::Vector< amrex::Geometry > &geom, const int current_N_level)
Function to calculate the relative B field error used in the predictor corrector loop.
Definition: Fields.cpp:1199
bool m_open_boundary
Definition: Fields.H:519
void FillBoundary(const amrex::Periodicity &period, const int lev, const int islice, Args...comps)
call amrex FillBoundary for multiple fields
Definition: Fields.H:471
void InSituComputeDiags(int step, amrex::Real time, int islice, const amrex::Geometry &geom3D, int max_step, amrex::Real max_time)
Definition: Fields.cpp:1259
void InitializeSlices(int lev, int islice, const amrex::Vector< amrex::Geometry > &geom)
Initialize all required fields to zero and interpolate from lev-1 to lev if needed.
Definition: Fields.cpp:545
void ShiftSlices(int lev)
Shift slices by 1 element: slices (1,2) are then stored in (2,3).
Definition: Fields.cpp:596
void AllocData(int lev, amrex::Geometry const &geom, const amrex::BoxArray &slice_ba, const amrex::DistributionMapping &slice_dm, int bin_size)
Definition: Fields.cpp:38
void setVal(const amrex::Real val, const int lev, const int islice, Args...comps)
set all selected fields to a value
Definition: Fields.H:309
amrex::Vector< amrex::MultiFab > m_slices
Definition: Fields.H:505
void LevelUpBoundary(amrex::Vector< amrex::Geometry > const &geom, const int lev, const int which_slice, const std::string &component, const amrex::IntVect outer_edge, const amrex::IntVect inner_edge)
Interpolate values from coarse grid (lev-1) to the boundary of the fine grid (lev)....
Definition: Fields.cpp:778
amrex::Vector< amrex::Real > m_insitu_rdata
Definition: Fields.H:530
bool m_do_symmetrize
Definition: Fields.H:501
void InitialBfieldGuess(const amrex::Real relative_Bfield_error, const amrex::Real predcorr_B_error_tolerance, const int lev)
Sets the initial guess of the B field from the two previous slices.
Definition: Fields.cpp:1120
amrex::Gpu::PinnedVector< amrex::Real > m_rel_z_vec_cpu
Definition: Fields.H:513
void shift(const int lev, const int islice_dst, const int islice_src, Args...comps)
copy all selected fields between slices
Definition: Fields.H:367
void InSituWriteToFile(int step, amrex::Real time, const amrex::Geometry &geom3D, int max_step, amrex::Real max_time)
Definition: Fields.cpp:1322
void add(const int lev, const int islice_dst, const char *const (&comps_dst)[ncomps], const int islice_src, const char *const (&comps_src)[ncomps])
add all selected fields between slices or on the same slice Uses references to C-arrays as argument s...
Definition: Fields.H:436
std::string m_insitu_file_prefix
Definition: Fields.H:534
void mult(const amrex::Real val, const int lev, const int islice, Args...comps)
multiply all selected fields with a value
Definition: Fields.H:338
void SetBoundaryCondition(amrex::Vector< amrex::Geometry > const &geom, const int lev, const int which_slice, std::string component, amrex::MultiFab &&staging_area)
Set up boundary conditions before poisson solve lev==0: leave at zero or add open boundaries lev>0: i...
Definition: Fields.cpp:685
amrex::Vector< std::unique_ptr< FFTPoissonSolver > > m_poisson_solver
Definition: Fields.H:103
Fields(const int nlev)
Definition: Fields.cpp:25
const amrex::MultiFab & getSlices(int lev) const
Definition: Fields.H:113
void SymmetrizeFields(int field_comp, const int lev, const int symm_x, const int symm_y)
Symmetrize fields by averaging over (x,y), symm_x*(-x,y), symm_y*(x,-y) and symm_x*symm_y*(-x,...
Definition: Fields.cpp:1084
static amrex::IntVect m_poisson_nguards
Definition: Fields.H:497
amrex::MultiFab getStagingArea(const int lev)
Definition: Fields.H:125
bool m_any_neutral_background
Definition: Fields.H:523
bool m_explicit
Definition: Fields.H:521
amrex::Vector< amrex::FArrayBox > & getTmpDensities()
Definition: Fields.H:129
amrex::MultiFab getField(const int lev, const int islice, const std::string comp)
Definition: Fields.H:119
void SolvePoissonPsiExmByEypBxEzBz(amrex::Vector< amrex::Geometry > const &geom, const int current_N_level)
Compute Psi, ExmBy, EypBx, Ez and Bz on the slice container from J by solving three poisson equations...
Definition: Fields.cpp:856
amrex::IntVect m_source_nguard
Definition: Fields.H:517
static constexpr int m_insitu_nrp
Definition: Fields.H:525
amrex::Gpu::DeviceVector< amrex::Real > m_rel_z_vec
Definition: Fields.H:511
void LevelUp(amrex::Vector< amrex::Geometry > const &geom, const int lev, const int which_slice, const std::string &component)
Interpolate the full field from the coarse grid (lev-1) to the fine grid (lev).
Definition: Fields.cpp:823
amrex::MultiFab & getSlices(int lev)
Definition: Fields.H:109
void AddRhoIons(const int lev)
Definition: Fields.cpp:614
bool m_do_dirichlet_poisson
Definition: Fields.H:507
bool m_extended_solve
Definition: Fields.H:499
amrex::Vector< amrex::Real > m_insitu_sum_rdata
Definition: Fields.H:532
amrex::Vector< amrex::FArrayBox > m_tmp_densities
Definition: Fields.H:509
void SolvePoissonEz(amrex::Vector< amrex::Geometry > const &geom, const int current_N_level, const int which_slice=WhichSlice::This)
Compute Ez on the slice container from J by solving a Poisson equation. This function does all the ne...
Definition: Fields.cpp:966
static amrex::IntVect m_slices_nguards
Definition: Fields.H:495
void MixAndShiftBfields(const amrex::Real relative_Bfield_error, const amrex::Real relative_Bfield_error_prev_iter, const amrex::Real predcorr_B_mixing_factor, const int lev)
Mixes the B field with the calculated current and previous iteration of it and shifts the current to ...
Definition: Fields.cpp:1144
void checkInit()
Definition: Fields.H:131
int m_insitu_period
Definition: Fields.H:528
amrex::Vector< amrex::MultiFab > & getSlices()
Definition: Fields.H:105
amrex::IntVect m_exmby_eypbx_nguard
Definition: Fields.H:515
void duplicate(const int lev, const int islice_dst, const char *const (&comps_dst)[ncomps], const int islice_src, const char *const (&comps_src)[ncomps])
copy all selected fields between slices or on the same slice Uses references to C-arrays as argument ...
Definition: Fields.H:399
void Copy(const int lev, const int i_slice, FieldDiagnosticData &fd, const amrex::Geometry &calc_geom, MultiLaser &multi_laser)
Copy between the full FArrayBox and slice MultiFab.
Definition: Fields.cpp:425
Singleton class, that intialize, runs and finalizes the simulation.
Definition: Hipace.H:70
Definition: MultiLaser.H:99
AMREX_GPU_HOST_DEVICE const IntVect & smallEnd() const &noexcept
AMREX_GPU_HOST_DEVICE const IntVect & bigEnd() const &noexcept
const Real * CellSize() const noexcept
void FillBoundary(bool cross=false)
Array4< typename FabArray< FArrayBox >::value_type const > array(const MFIter &mfi) const noexcept
const Real * ProbLo() const noexcept
const Real * ProbHi() const noexcept
int i
Definition: MakeOpenBoundary.py:152
j
Definition: MakeOpenBoundary.py:128
AMREX_GPU_HOST_DEVICE Long at(T const &, Long offset) noexcept
std::ostream & ErrorStream()
std::enable_if_t< std::is_integral< T >::value > ParallelFor(TypeList< CTOs... >, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
const int[]
str
Definition: checksumAPI.py:112
Definition: GPUUtil.H:98
Direction of each dimension. Can be used for clean handling 2D vs. 3D in the future.
Definition: Fields.H:59
dir
Definition: Fields.H:60
@ z
Definition: Fields.H:60
@ x
Definition: Fields.H:60
@ y
Definition: Fields.H:60
This struct holds data for one field diagnostic on one MR level.
Definition: Diagnostic.H:21
describes which slice with respect to the currently calculated is used
Definition: Fields.H:27
slice
Definition: Fields.H:28
@ RhomJzIons
Definition: Fields.H:28
@ PCPrevIter
Definition: Fields.H:28
@ Next
Definition: Fields.H:28
@ Salame
Definition: Fields.H:28
@ Previous
Definition: Fields.H:28
@ This
Definition: Fields.H:28
@ PCIter
Definition: Fields.H:28
@ N
Definition: Fields.H:28
Definition: Fields.H:31
void multi_emplace(int &n, Args...comps)
Definition: Fields.H:45
int operator[](const std::string &str) const
Definition: Fields.H:32