WarpX
VectorPoissonSolver.H
Go to the documentation of this file.
1 /* Copyright 2022 S. Eric Clark, LLNL
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef ABLASTR_VECTOR_POISSON_SOLVER_H
8 #define ABLASTR_VECTOR_POISSON_SOLVER_H
9 
10 #include <ablastr/constant.H>
13 #include <ablastr/utils/TextMsg.H>
15 
16 #include <AMReX_Array.H>
17 #include <AMReX_Array4.H>
18 #include <AMReX_BLassert.H>
19 #include <AMReX_Box.H>
20 #include <AMReX_BoxArray.H>
21 #include <AMReX_Config.H>
23 #include <AMReX_FArrayBox.H>
24 #include <AMReX_FabArray.H>
25 #include <AMReX_Geometry.H>
26 #include <AMReX_GpuControl.H>
27 #include <AMReX_GpuLaunch.H>
28 #include <AMReX_GpuQualifiers.H>
29 #include <AMReX_IndexType.H>
30 #include <AMReX_IntVect.H>
31 #include <AMReX_LO_BCTYPES.H>
32 #include <AMReX_MFIter.H>
33 #include <AMReX_MFInterp_C.H>
35 #include <AMReX_MLLinOp.H>
36 #include <AMReX_MLMG.H>
37 #include <AMReX_MultiFab.H>
38 #include <AMReX_Parser.H>
39 #include <AMReX_REAL.H>
40 #include <AMReX_SPACE.H>
41 #include <AMReX_Vector.H>
42 #ifdef AMREX_USE_EB
43 # include <AMReX_EBFabFactory.H>
44 #endif
45 
46 #include <array>
47 #include <optional>
48 
49 namespace ablastr::fields {
50 
80 template<
81  typename T_BoundaryHandler,
82  typename T_PostACalculationFunctor = std::nullopt_t,
83  typename T_FArrayBoxFactory = void
84 >
85 void
88  amrex::Real const relative_tolerance,
89  amrex::Real absolute_tolerance,
90  int const max_iters,
91  int const verbosity,
94  amrex::Vector<amrex::BoxArray> const& grids,
95  T_BoundaryHandler const boundary_handler,
96  bool const do_single_precision_comms = false,
97  std::optional<amrex::Vector<amrex::IntVect> > rel_ref_ratio = std::nullopt,
98  [[maybe_unused]] T_PostACalculationFunctor post_A_calculation = std::nullopt,
99  [[maybe_unused]] std::optional<amrex::Real const> current_time = std::nullopt, // only used for EB
100  [[maybe_unused]] std::optional<amrex::Vector<T_FArrayBoxFactory const *> > eb_farray_box_factory = std::nullopt // only used for EB
101 )
102 {
103  using namespace amrex::literals;
104 
105  if (!rel_ref_ratio.has_value()) {
106  ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE(curr.size() == 1u,
107  "rel_ref_ratio must be set if mesh-refinement is used");
108  rel_ref_ratio = amrex::Vector<amrex::IntVect>{{amrex::IntVect(AMREX_D_DECL(1, 1, 1))}};
109  }
110 
111  auto const finest_level = static_cast<int>(curr.size()) - 1;
112 
113  // scale J appropriately; also determine if current is zero everywhere
114  amrex::Real max_comp_J = 0.0;
115  for (int lev=0; lev<=finest_level; lev++) {
116  for (int adim=0; adim<3; adim++) {
117  curr[lev][adim]->mult(-1._rt*ablastr::constant::SI::mu0); // Unscaled below
118  max_comp_J = amrex::max(max_comp_J, curr[lev][adim]->norm0());
119  }
120  }
122 
123  const bool always_use_bnorm = (max_comp_J > 0);
124  if (!always_use_bnorm) {
125  if (absolute_tolerance == 0.0) { absolute_tolerance = amrex::Real(1e-6); }
127  "MagnetostaticSolver",
128  "Max norm of J is 0",
130  );
131  }
132 
133  const amrex::LPInfo& info = amrex::LPInfo();
134 
135  // Loop over dimensions of A to solve each component individually
136  for (int lev=0; lev<=finest_level; lev++) {
138  {geom[lev]}, {grids[lev]}, {dmap[lev]}, info
139 #if defined(AMREX_USE_EB)
140  , {eb_farray_box_factory.value()[lev]}
141 #endif
142  );
144  {geom[lev]}, {grids[lev]}, {dmap[lev]}, info
145 #if defined(AMREX_USE_EB)
146  , {eb_farray_box_factory.value()[lev]}
147 #endif
148  );
150  {geom[lev]}, {grids[lev]}, {dmap[lev]}, info
151 #if defined(AMREX_USE_EB)
152  , {eb_farray_box_factory.value()[lev]}
153 #endif
154  );
155 
156  amrex::Array<amrex::MLEBNodeFDLaplacian*,3> linop = {&linopx,&linopy,&linopz};
158 
159  for (int adim=0; adim<3; adim++) {
160  // Solve the Poisson equation
161  // This is solving the self fields using the magnetostatic solver in the lab frame
162 
163  // Note: this assumes that beta is zero
164  linop[adim]->setSigma({AMREX_D_DECL(1._rt, 1._rt, 1._rt)});
165 
166 #if defined(AMREX_USE_EB)
167  // Set Homogeneous Dirichlet Boundary on EB
168  linop[adim]->setEBDirichlet(0_rt);
169 #endif
170 
171 #ifdef WARPX_DIM_RZ
172  linop[adim]->setRZ(true);
173 #endif
174 
175  linop[adim]->setDomainBC( boundary_handler.lobc[adim], boundary_handler.hibc[adim] );
176 
177  mlmg[adim] = std::make_unique<amrex::MLMG>(*linop[adim]);
178 
179  mlmg[adim]->setVerbose(verbosity);
180  mlmg[adim]->setMaxIter(max_iters);
181  mlmg[adim]->setAlwaysUseBNorm(always_use_bnorm);
182 
183  // Solve Poisson equation at lev
184  mlmg[adim]->solve( {A[lev][adim]}, {curr[lev][adim]},
185  relative_tolerance, absolute_tolerance );
186 
187  // Synchronize the ghost cells, do halo exchange
189  *A[lev][adim], A[lev][adim]->nGrowVect(),
190  do_single_precision_comms,
191  geom[lev].periodicity(), false);
192 
193  // needed for solving the levels by levels:
194  // - coarser level is initial guess for finer level
195  // - coarser level provides boundary values for finer level patch
196  // Interpolation from phi[lev] to phi[lev+1]
197  // (This provides both the boundary conditions and initial guess for phi[lev+1])
198  if (lev < finest_level) {
199 
200  // Allocate A_cp for lev+1
201  amrex::BoxArray ba = A[lev+1][adim]->boxArray();
202  const amrex::IntVect& refratio = rel_ref_ratio.value()[lev];
203  ba.coarsen(refratio);
204  const int ncomp = linop[adim]->getNComp();
205  amrex::MultiFab A_cp(ba, A[lev+1][adim]->DistributionMap(), ncomp, 1);
206 
207  // Copy from A[lev] to A_cp (in parallel)
209  const amrex::Periodicity& crse_period = geom[lev].periodicity();
210 
212  A_cp,
213  *A[lev][adim],
214  0,
215  0,
216  1,
217  ng,
218  ng,
219  do_single_precision_comms,
220  crse_period
221  );
222 
223  // Local interpolation from A_cp to A[lev+1]
224 #ifdef AMREX_USE_OMP
225 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
226 #endif
227  for (amrex::MFIter mfi(*A[lev+1][adim],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
228  {
229  amrex::Array4<amrex::Real> const& A_fp_arr = A[lev+1][adim]->array(mfi);
230  amrex::Array4<amrex::Real> const& A_cp_arr = A_cp.array(mfi);
231 
232  details::PoissonInterpCPtoFP const interp(A_fp_arr, A_cp_arr, refratio);
233 
234  amrex::Box const b = mfi.tilebox(A[lev + 1][adim]->ixType().toIntVect());
236  }
237 
238  }
239 
240  // Unscale current
241  curr[lev][adim]->mult(-1._rt/ablastr::constant::SI::mu0);
242  } // Loop over adim
243  // Run additional operations, such as calculation of the B fields for embedded boundaries
245  if (post_A_calculation.has_value()) {
246  post_A_calculation.value()(mlmg, lev);
247  }
248  }
249  } // loop over lev(els)
250 }
251 } // namepace Magnetostatic
252 #endif
#define AMREX_D_DECL(a, b, c)
#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
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
static constexpr auto mu0
vacuum permeability: magnetic permeability of vacuum = 4.0e-7 * pi [H/m]
Definition: constant.H:48
Definition: IntegratedGreenFunctionSolver.cpp:33
void computeVectorPotential(amrex::Vector< amrex::Array< amrex::MultiFab *, 3 > > const &curr, amrex::Vector< amrex::Array< amrex::MultiFab *, 3 > > &A, 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, bool const do_single_precision_comms=false, std::optional< amrex::Vector< amrex::IntVect > > rel_ref_ratio=std::nullopt, [[maybe_unused]] T_PostACalculationFunctor post_A_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: VectorPoissonSolver.H:86
void FillBoundary(amrex::MultiFab &mf, amrex::IntVect ng, bool do_single_precision_comms, const amrex::Periodicity &period, std::optional< bool > nodal_sync)
Definition: Communication.cpp:71
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)
IntVect nGrowVect(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
value
Definition: updateAMReX.py:141
int verbosity()