7 #ifndef ABLASTR_POISSON_SOLVER_H
8 #define ABLASTR_POISSON_SOLVER_H
18 #if defined(WARPX_USE_PSATD) && defined(WARPX_DIM_3D)
27 #include <AMReX_Config.H>
48 #if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)
93 typename T_BoundaryHandler,
94 typename T_PostPhiCalculationFunctor = std::nullopt_t,
95 typename T_FArrayBoxFactory =
void
100 std::array<amrex::Real, 3>
const beta,
101 amrex::Real
const relative_tolerance,
102 amrex::Real absolute_tolerance,
108 T_BoundaryHandler
const boundary_handler,
109 [[maybe_unused]]
bool is_solver_multigrid,
110 bool const do_single_precision_comms =
false,
112 [[maybe_unused]] T_PostPhiCalculationFunctor post_phi_calculation = std::nullopt,
113 [[maybe_unused]] std::optional<amrex::Real const> current_time = std::nullopt,
117 using namespace amrex::literals;
121 if (!rel_ref_ratio.has_value()) {
123 "rel_ref_ratio must be set if mesh-refinement is used");
127 auto const finest_level =
static_cast<int>(rho.
size() - 1);
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());
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",
146 #if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ))
152 for (
int lev=0; lev<=finest_level; lev++) {
155 #if defined(WARPX_DIM_1D_Z)
157 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
163 #if (defined(WARPX_USE_PSATD) && defined(WARPX_DIM_3D))
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 ) {
182 rho[lev]->mult(-1._rt/
ep0);
184 #if !(defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ))
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]));
201 if (max_semicoarsening_level > 0) {
208 #if defined(AMREX_USE_EB) || defined(WARPX_DIM_RZ)
213 #if defined(AMREX_USE_EB)
214 , {eb_farray_box_factory.value()[lev]}
221 #if defined(WARPX_DIM_RZ)
222 linop.
setSigma({0._rt, 1._rt-beta_solver[1]*beta_solver[1]});
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])});
230 #if defined(AMREX_USE_EB)
233 if (boundary_handler.phi_EB_only_t) {
234 linop.setEBDirichlet(boundary_handler.potential_eb_t(current_time.value()));
237 linop.setEBDirichlet(boundary_handler.getPhiEB(current_time.value()));
248 linop.setDomainBC( boundary_handler.lobc, boundary_handler.hibc );
258 mlmg.
solve( {phi[lev]}, {rho[lev]},
259 relative_tolerance, absolute_tolerance );
266 if (lev < finest_level) {
272 const int ncomp = linop.getNComp();
287 do_single_precision_comms,
293 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
301 amrex::Box const b = mfi.tilebox(phi[lev + 1]->ixType().toIntVect());
309 if (post_phi_calculation.has_value()) {
310 post_phi_calculation.value()(mlmg, lev);
#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)
beta
Definition: stencil.py:434
value
Definition: updateAMReX.py:141
Definition: Interpolate.H:33
LPInfo & setSemicoarsening(bool x) noexcept
LPInfo & setMaxSemicoarseningLevel(int n) noexcept
LPInfo & setSemicoarseningDirection(int n) noexcept