WarpX
PoissonSolver.H
Go to the documentation of this file.
1 /* Copyright 2019-2022 Axel Huebl, Remi Lehe
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef ABLASTR_POISSON_SOLVER_H
8 #define ABLASTR_POISSON_SOLVER_H
9 
10 #include <ablastr/constant.H>
12 #include <ablastr/utils/TextMsg.H>
17 
18 #if defined(WARPX_USE_PSATD) && defined(WARPX_DIM_3D)
20 #endif
21 
22 #include <AMReX_Array.H>
23 #include <AMReX_Array4.H>
24 #include <AMReX_BLassert.H>
25 #include <AMReX_Box.H>
26 #include <AMReX_BoxArray.H>
27 #include <AMReX_Config.H>
29 #include <AMReX_FArrayBox.H>
30 #include <AMReX_FabArray.H>
31 #include <AMReX_Geometry.H>
32 #include <AMReX_GpuControl.H>
33 #include <AMReX_GpuLaunch.H>
34 #include <AMReX_GpuQualifiers.H>
35 #include <AMReX_IndexType.H>
36 #include <AMReX_IntVect.H>
37 #include <AMReX_LO_BCTYPES.H>
38 #include <AMReX_MFIter.H>
39 #include <AMReX_MFInterp_C.H>
40 #include <AMReX_MLMG.H>
41 #include <AMReX_MLLinOp.H>
43 #include <AMReX_MultiFab.H>
44 #include <AMReX_Parser.H>
45 #include <AMReX_REAL.H>
46 #include <AMReX_SPACE.H>
47 #include <AMReX_Vector.H>
48 #if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)
50 #endif
51 #ifdef AMREX_USE_EB
52 # include <AMReX_EBFabFactory.H>
53 #endif
54 
55 #include <array>
56 #include <optional>
57 
58 
59 namespace ablastr::fields {
60 
92 template<
93  typename T_BoundaryHandler,
94  typename T_PostPhiCalculationFunctor = std::nullopt_t,
95  typename T_FArrayBoxFactory = void
96 >
97 void
100  std::array<amrex::Real, 3> const beta,
101  amrex::Real const relative_tolerance,
102  amrex::Real absolute_tolerance,
103  int const max_iters,
104  int const verbosity,
105  amrex::Vector<amrex::Geometry> const& geom,
107  amrex::Vector<amrex::BoxArray> const& grids,
108  T_BoundaryHandler const boundary_handler,
109  [[maybe_unused]] bool is_solver_multigrid,
110  bool const do_single_precision_comms = false,
111  std::optional<amrex::Vector<amrex::IntVect> > rel_ref_ratio = std::nullopt,
112  [[maybe_unused]] T_PostPhiCalculationFunctor post_phi_calculation = std::nullopt,
113  [[maybe_unused]] std::optional<amrex::Real const> current_time = std::nullopt, // only used for EB
114  [[maybe_unused]] std::optional<amrex::Vector<T_FArrayBoxFactory const *> > eb_farray_box_factory = std::nullopt // only used for EB
115 )
116 {
117  using namespace amrex::literals;
118 
119  ABLASTR_PROFILE("computePhi", false);
120 
121  if (!rel_ref_ratio.has_value()) {
123  "rel_ref_ratio must be set if mesh-refinement is used");
124  rel_ref_ratio = amrex::Vector<amrex::IntVect>{{amrex::IntVect(AMREX_D_DECL(1, 1, 1))}};
125  }
126 
127  auto const finest_level = static_cast<int>(rho.size() - 1);
128 
129  // determine if rho is zero everywhere
130  amrex::Real max_norm_b = 0.0;
131  for (int lev=0; lev<=finest_level; lev++) {
132  max_norm_b = amrex::max(max_norm_b, rho[lev]->norm0());
133  }
135 
136  const bool always_use_bnorm = (max_norm_b > 0);
137  if (!always_use_bnorm) {
138  if (absolute_tolerance == 0.0) { absolute_tolerance = amrex::Real(1e-6); }
140  "ElectrostaticSolver",
141  "Max norm of rho is 0",
143  );
144  }
145 
146 #if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ))
147  amrex::LPInfo info;
148 #else
149  const amrex::LPInfo info;
150 #endif
151 
152  for (int lev=0; lev<=finest_level; lev++) {
153  // Set the value of beta
155 #if defined(WARPX_DIM_1D_Z)
156  {{ beta[2] }}; // beta_x and beta_z
157 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
158  {{ beta[0], beta[2] }}; // beta_x and beta_z
159 #else
160  {{ beta[0], beta[1], beta[2] }};
161 #endif
162 
163 #if (defined(WARPX_USE_PSATD) && defined(WARPX_DIM_3D))
164  // Use the Integrated Green Function solver (FFT) on the coarsest level if it was selected
165  if(!is_solver_multigrid && lev==0){
167  {AMREX_D_DECL(geom[lev].CellSize(0)/std::sqrt(1._rt-beta_solver[0]*beta_solver[0]),
168  geom[lev].CellSize(1)/std::sqrt(1._rt-beta_solver[1]*beta_solver[1]),
169  geom[lev].CellSize(2)/std::sqrt(1._rt-beta_solver[2]*beta_solver[2]))};
170  if ( max_norm_b == 0 ) {
171  phi[lev]->setVal(0);
172  } else {
173  computePhiIGF( *rho[lev], *phi[lev], dx_igf, grids[lev] );
174  }
175  continue;
176  }
177 #endif
178 
179  // Use the Multigrid (MLMG) solver if selected or on refined patches
180  // but first scale rho appropriately
181  using namespace ablastr::constant::SI;
182  rho[lev]->mult(-1._rt/ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect!
183 
184 #if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ))
185  // Determine whether to use semi-coarsening
187  {AMREX_D_DECL(geom[lev].CellSize(0)/std::sqrt(1._rt-beta_solver[0]*beta_solver[0]),
188  geom[lev].CellSize(1)/std::sqrt(1._rt-beta_solver[1]*beta_solver[1]),
189  geom[lev].CellSize(2)/std::sqrt(1._rt-beta_solver[2]*beta_solver[2]))};
190  int max_semicoarsening_level = 0;
191  int semicoarsening_direction = -1;
192  const auto min_dir = static_cast<int>(std::distance(dx_scaled.begin(),
193  std::min_element(dx_scaled.begin(),dx_scaled.end())));
194  const auto max_dir = static_cast<int>(std::distance(dx_scaled.begin(),
195  std::max_element(dx_scaled.begin(),dx_scaled.end())));
196  if (dx_scaled[max_dir] > dx_scaled[min_dir]) {
197  semicoarsening_direction = max_dir;
198  max_semicoarsening_level = static_cast<int>
199  (std::log2(dx_scaled[max_dir]/dx_scaled[min_dir]));
200  }
201  if (max_semicoarsening_level > 0) {
202  info.setSemicoarsening(true);
203  info.setMaxSemicoarseningLevel(max_semicoarsening_level);
204  info.setSemicoarseningDirection(semicoarsening_direction);
205  }
206 #endif
207 
208 #if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)
209  // In the presence of EB or RZ: the solver assumes that the beam is
210  // propagating along one of the axes of the grid, i.e. that only *one*
211  // of the components of `beta` is non-negligible.
212  amrex::MLEBNodeFDLaplacian linop( {geom[lev]}, {grids[lev]}, {dmap[lev]}, info
213 #if defined(AMREX_USE_EB)
214  , {eb_farray_box_factory.value()[lev]}
215 #endif
216  );
217 
218  // Note: this assumes that the beam is propagating along
219  // one of the axes of the grid, i.e. that only *one* of the
220  // components of `beta` is non-negligible. // we use this
221 #if defined(WARPX_DIM_RZ)
222  linop.setSigma({0._rt, 1._rt-beta_solver[1]*beta_solver[1]});
223 #else
224  linop.setSigma({AMREX_D_DECL(
225  1._rt-beta_solver[0]*beta_solver[0],
226  1._rt-beta_solver[1]*beta_solver[1],
227  1._rt-beta_solver[2]*beta_solver[2])});
228 #endif
229 
230 #if defined(AMREX_USE_EB)
231  // if the EB potential only depends on time, the potential can be passed
232  // as a float instead of a callable
233  if (boundary_handler.phi_EB_only_t) {
234  linop.setEBDirichlet(boundary_handler.potential_eb_t(current_time.value()));
235  }
236  else
237  linop.setEBDirichlet(boundary_handler.getPhiEB(current_time.value()));
238 #endif
239 #else
240  // In the absence of EB and RZ: use a more generic solver
241  // that can handle beams propagating in any direction
242  amrex::MLNodeTensorLaplacian linop( {geom[lev]}, {grids[lev]},
243  {dmap[lev]}, info );
244  linop.setBeta( beta_solver ); // for the non-axis-aligned solver
245 #endif
246 
247  // Solve the Poisson equation
248  linop.setDomainBC( boundary_handler.lobc, boundary_handler.hibc );
249 #ifdef WARPX_DIM_RZ
250  linop.setRZ(true);
251 #endif
252  amrex::MLMG mlmg(linop); // actual solver defined here
253  mlmg.setVerbose(verbosity);
254  mlmg.setMaxIter(max_iters);
255  mlmg.setAlwaysUseBNorm(always_use_bnorm);
256 
257  // Solve Poisson equation at lev
258  mlmg.solve( {phi[lev]}, {rho[lev]},
259  relative_tolerance, absolute_tolerance );
260 
261  // needed for solving the levels by levels:
262  // - coarser level is initial guess for finer level
263  // - coarser level provides boundary values for finer level patch
264  // Interpolation from phi[lev] to phi[lev+1]
265  // (This provides both the boundary conditions and initial guess for phi[lev+1])
266  if (lev < finest_level) {
267 
268  // Allocate phi_cp for lev+1
269  amrex::BoxArray ba = phi[lev+1]->boxArray();
270  const amrex::IntVect& refratio = rel_ref_ratio.value()[lev];
271  ba.coarsen(refratio);
272  const int ncomp = linop.getNComp();
273  amrex::MultiFab phi_cp(ba, phi[lev+1]->DistributionMap(), ncomp, 1);
274 
275  // Copy from phi[lev] to phi_cp (in parallel)
277  const amrex::Periodicity& crse_period = geom[lev].periodicity();
278 
280  phi_cp,
281  *phi[lev],
282  0,
283  0,
284  1,
285  ng,
286  ng,
287  do_single_precision_comms,
288  crse_period
289  );
290 
291  // Local interpolation from phi_cp to phi[lev+1]
292 #ifdef AMREX_USE_OMP
293 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
294 #endif
295  for (amrex::MFIter mfi(*phi[lev + 1], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
296  amrex::Array4<amrex::Real> const phi_fp_arr = phi[lev + 1]->array(mfi);
297  amrex::Array4<amrex::Real const> const phi_cp_arr = phi_cp.array(mfi);
298 
299  details::PoissonInterpCPtoFP const interp(phi_fp_arr, phi_cp_arr, refratio);
300 
301  amrex::Box const b = mfi.tilebox(phi[lev + 1]->ixType().toIntVect());
303  }
304 
305  }
306 
307  // Run additional operations, such as calculation of the E field for embedded boundaries
309  if (post_phi_calculation.has_value()) {
310  post_phi_calculation.value()(mlmg, lev);
311  }
312  }
313 
314  } // loop over lev(els)
315 }
316 
317 } // namespace ablastr::fields
318 
319 #endif // ABLASTR_POISSON_SOLVER_H
#define AMREX_D_DECL(a, b, c)
#define ABLASTR_PROFILE(fname, sync)
Definition: ProfilerWrapper.H:55
#define ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition: TextMsg.H:75
BoxArray & coarsen(int refinement_ratio)
Array4< typename FabArray< FArrayBox >::value_type const > array(const MFIter &mfi) const noexcept
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE IntVect TheUnitVector() noexcept
void setSigma(Array< Real, AMREX_SPACEDIM > const &a_sigma) noexcept
void setVerbose(int v) noexcept
void setAlwaysUseBNorm(int flag) noexcept
RT solve(const Vector< AMF * > &a_sol, const Vector< AMF const * > &a_rhs, RT a_tol_rel, RT a_tol_abs, const char *checkpoint_file=nullptr)
void setMaxIter(int n) noexcept
void setBeta(Array< Real, AMREX_SPACEDIM > const &a_beta) noexcept
Long size() const noexcept
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void interp(int j, int k, int l, amrex::Array4< amrex::Real > const &fine, amrex::Array4< amrex::Real const > const &crse, const amrex::IntVect r_ratio, IntVect const &type) noexcept
Definition: Interpolate_K.H:9
Definition: constant.H:40
static constexpr auto ep0
vacuum permittivity: dielectric permittivity of vacuum [F/m]
Definition: constant.H:46
Definition: IntegratedGreenFunctionSolver.cpp:33
void computePhiIGF(amrex::MultiFab const &rho, amrex::MultiFab &phi, std::array< amrex::Real, 3 > const &cell_size, amrex::BoxArray const &ba)
Compute the electrostatic potential using the Integrated Green Function method as in http://dx....
Definition: IntegratedGreenFunctionSolver.cpp:36
void computePhi(amrex::Vector< amrex::MultiFab * > const &rho, amrex::Vector< amrex::MultiFab * > &phi, std::array< amrex::Real, 3 > const beta, amrex::Real const relative_tolerance, amrex::Real absolute_tolerance, int const max_iters, int const verbosity, amrex::Vector< amrex::Geometry > const &geom, amrex::Vector< amrex::DistributionMapping > const &dmap, amrex::Vector< amrex::BoxArray > const &grids, T_BoundaryHandler const boundary_handler, [[maybe_unused]] bool is_solver_multigrid, bool const do_single_precision_comms=false, std::optional< amrex::Vector< amrex::IntVect > > rel_ref_ratio=std::nullopt, [[maybe_unused]] T_PostPhiCalculationFunctor post_phi_calculation=std::nullopt, [[maybe_unused]] std::optional< amrex::Real const > current_time=std::nullopt, [[maybe_unused]] std::optional< amrex::Vector< T_FArrayBoxFactory const * > > eb_farray_box_factory=std::nullopt)
Definition: PoissonSolver.H:98
void ParallelCopy(amrex::MultiFab &dst, const amrex::MultiFab &src, int src_comp, int dst_comp, int num_comp, const amrex::IntVect &src_nghost, const amrex::IntVect &dst_nghost, bool do_single_precision_comms, const amrex::Periodicity &period, amrex::FabArrayBase::CpOp op)
Definition: Communication.cpp:28
void WMRecordWarning(const std::string &topic, const std::string &text, const WarnPriority &priority=WarnPriority::medium)
Helper function to abbreviate the call to WarnManager::GetInstance().RecordWarning (recording a warni...
Definition: WarnManager.cpp:318
void ReduceRealMax(Vector< std::reference_wrapper< Real > > const &)
DistributionMapping const & DistributionMap(FabArrayBase const &fa)
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
bool TilingIfNotGPU() noexcept
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... >, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
std::array< T, N > Array
beta
Definition: stencil.py:434
value
Definition: updateAMReX.py:141
LPInfo & setSemicoarsening(bool x) noexcept
LPInfo & setMaxSemicoarseningLevel(int n) noexcept
LPInfo & setSemicoarseningDirection(int n) noexcept
int verbosity()