From d4c01c3dc434f5eea333e103a1d6a49670bae874 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Mon, 4 Nov 2024 16:37:21 -0800 Subject: [PATCH] FFT: OpenBC Solver This implements the algorithm of Hockney, Methods Comp. Phys. 9 (1970), 136-210 for solving Poisson's equation with open boundaries. --- Src/Base/AMReX_BoxArray.H | 3 +- Src/Base/AMReX_BoxArray.cpp | 21 ++++- Src/FFT/AMReX_FFT.H | 1 + Src/FFT/AMReX_FFT_OpenBCSolver.H | 142 +++++++++++++++++++++++++++++++ Src/FFT/AMReX_FFT_Poisson.H | 80 ++++++++++++++++- Src/FFT/AMReX_FFT_R2C.H | 43 ++++++---- Src/FFT/CMakeLists.txt | 1 + Src/FFT/Make.package | 3 +- Tests/FFT/OpenBC/CMakeLists.txt | 13 +++ Tests/FFT/OpenBC/GNUmakefile | 26 ++++++ Tests/FFT/OpenBC/Make.package | 1 + Tests/FFT/OpenBC/main.cpp | 125 +++++++++++++++++++++++++++ 12 files changed, 437 insertions(+), 22 deletions(-) create mode 100644 Src/FFT/AMReX_FFT_OpenBCSolver.H create mode 100644 Tests/FFT/OpenBC/CMakeLists.txt create mode 100644 Tests/FFT/OpenBC/GNUmakefile create mode 100644 Tests/FFT/OpenBC/Make.package create mode 100644 Tests/FFT/OpenBC/main.cpp diff --git a/Src/Base/AMReX_BoxArray.H b/Src/Base/AMReX_BoxArray.H index e85946872ce..066e3b073f7 100644 --- a/Src/Base/AMReX_BoxArray.H +++ b/Src/Base/AMReX_BoxArray.H @@ -69,7 +69,8 @@ namespace amrex */ [[nodiscard]] BoxArray decompose (Box const& domain, int nboxes, Array const& decomp - = {AMREX_D_DECL(true,true,true)}); + = {AMREX_D_DECL(true,true,true)}, + bool no_overlap = false); struct BARef { diff --git a/Src/Base/AMReX_BoxArray.cpp b/Src/Base/AMReX_BoxArray.cpp index 576d4cb870d..2e30a0a99c2 100644 --- a/Src/Base/AMReX_BoxArray.cpp +++ b/Src/Base/AMReX_BoxArray.cpp @@ -1891,7 +1891,7 @@ bool match (const BoxArray& x, const BoxArray& y) } BoxArray decompose (Box const& domain, int nboxes, - Array const& decomp) + Array const& decomp, bool no_overlap) { auto ndecomp = std::count(decomp.begin(), decomp.end(), true); @@ -2048,9 +2048,24 @@ BoxArray decompose (Box const& domain, int nboxes, ilo += domlo[0]; ihi += domlo[0]; Box b{IntVect(AMREX_D_DECL(ilo,jlo,klo)), - IntVect(AMREX_D_DECL(ihi,jhi,khi))}; + IntVect(AMREX_D_DECL(ihi,jhi,khi)), ixtyp}; if (b.ok()) { - bl.push_back(b.convert(ixtyp)); + if (no_overlap) { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (ixtyp.nodeCentered() && + b.bigEnd(idim) == ccdomain.bigEnd(idim)) + { + b.growHi(idim, 1); + } + } + } else { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (ixtyp.nodeCentered()) { + b.growHi(idim, 1); + } + } + } + bl.push_back(b); } AMREX_D_TERM(},},}) diff --git a/Src/FFT/AMReX_FFT.H b/Src/FFT/AMReX_FFT.H index 22bcb1413ec..cd7c0984c59 100644 --- a/Src/FFT/AMReX_FFT.H +++ b/Src/FFT/AMReX_FFT.H @@ -2,6 +2,7 @@ #define AMREX_FFT_H_ #include +#include #include #include diff --git a/Src/FFT/AMReX_FFT_OpenBCSolver.H b/Src/FFT/AMReX_FFT_OpenBCSolver.H new file mode 100644 index 00000000000..06e3806a5b6 --- /dev/null +++ b/Src/FFT/AMReX_FFT_OpenBCSolver.H @@ -0,0 +1,142 @@ +#ifndef AMREX_FFT_OPENBC_SOlVER_H_ +#define AMREX_FFT_OPENBC_SOlVER_H_ + +#include + +#include + +namespace amrex::FFT +{ + +template +class OpenBCSolver +{ +public: + using MF = typename R2C::MF; + using cMF = typename R2C::cMF; + + explicit OpenBCSolver (Box const& domain); + + template + void setGreensFunction (F const& greens_function); + + void solve (MF& phi, MF const& rho); + + [[nodiscard]] Box const& Domain () const { return m_domain; } + +private: + Box m_domain; + R2C m_r2c; + cMF m_G_fft; +}; + +template +OpenBCSolver::OpenBCSolver (Box const& domain) + : m_domain(domain), + m_r2c(Box(domain.smallEnd(), domain.bigEnd()+domain.length(), domain.ixType())) +{ + auto [sd, ord] = m_r2c.getSpectralData(); + amrex::ignore_unused(ord); + m_G_fft.define(sd->boxArray(), sd->DistributionMap(), 1, 0); +} + +template +template +void OpenBCSolver::setGreensFunction (F const& greens_function) +{ + auto* infab = detail::get_fab(m_r2c.m_rx); + auto const& lo = m_domain.smallEnd(); + auto const& lo3 = lo.dim3(); + auto const& len = m_domain.length3d(); + if (infab) { + auto const& a = infab->array(); + auto box = infab->box(); + GpuArray nimages{1,1,1}; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (box.smallEnd(idim) == lo[idim] && box.length(idim) == 2*len[idim]) { + box.growHi(idim, -len[idim]+1); // +1 to include the middle plane + nimages[idim] = 2; + } + } + AMREX_ASSERT(nimages[0] == 2); + box.shift(-lo); + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (i == len[0] || j == len[1] || k == len[2]) { + a(i+lo3.x,j+lo3.y,k+lo3.z) = T(0); + } else { + auto ii = i; + auto jj = (j > len[1]) ? 2*len[1]-j : j; + auto kk = (k > len[2]) ? 2*len[2]-k : k; + auto G = greens_function(ii+lo3.x,jj+lo3.y,kk+lo3.z); + for (int koff = 0; koff < nimages[2]; ++koff) { + int k2 = (koff == 0) ? k : 2*len[2]-k; + if (k2 == 2*len[2]) { continue; } + for (int joff = 0; joff < nimages[1]; ++joff) { + int j2 = (joff == 0) ? j : 2*len[1]-j; + if (j2 == 2*len[1]) { continue; } + for (int ioff = 0; ioff < nimages[0]; ++ioff) { + int i2 = (ioff == 0) ? i : 2*len[0]-i; + if (i2 == 2*len[0]) { continue; } + a(i2+lo3.x,j2+lo3.y,k2+lo3.z) = G; + } + } + } + } + }); + } + + m_r2c.forward(m_r2c.m_rx); + + auto [sd, ord] = m_r2c.getSpectralData(); + amrex::ignore_unused(ord); + auto const* srcfab = detail::get_fab(*sd); + if (srcfab) { + auto* dstfab = detail::get_fab(m_G_fft); + if (dstfab) { +#if defined(AMREX_USE_GPU) + Gpu::dtod_memcpy_async +#else + std::memcpy +#endif + (dstfab->dataPtr(), srcfab->dataPtr(), dstfab->nBytes()); + } else { + amrex::Abort("FFT::OpenBCSolver: how did this happen"); + } + } +} + +template +void OpenBCSolver::solve (MF& phi, MF const& rho) +{ + auto& inmf = m_r2c.m_rx; + inmf.setVal(T(0)); + inmf.ParallelCopy(rho, 0, 0, 1); + + m_r2c.forward(inmf); + + auto scaling_factor = T(1) / T(m_r2c.m_real_domain.numPts()); + + auto const* gfab = detail::get_fab(m_G_fft); + if (gfab) { + auto [sd, ord] = m_r2c.getSpectralData(); + amrex::ignore_unused(ord); + auto* rhofab = detail::get_fab(*sd); + if (rhofab) { + auto* pdst = rhofab->dataPtr(); + auto const* psrc = gfab->dataPtr(); + amrex::ParallelFor(rhofab->box().numPts(), [=] AMREX_GPU_DEVICE (Long i) + { + pdst[i] *= psrc[i] * scaling_factor; + }); + } else { + amrex::Abort("FFT::OpenBCSolver::solve: how did this happen?"); + } + } + + m_r2c.backward_doit(phi, phi.nGrowVect()); +} + +} + +#endif diff --git a/Src/FFT/AMReX_FFT_Poisson.H b/Src/FFT/AMReX_FFT_Poisson.H index b89ffcc85c0..062cd457d88 100644 --- a/Src/FFT/AMReX_FFT_Poisson.H +++ b/Src/FFT/AMReX_FFT_Poisson.H @@ -8,7 +8,8 @@ namespace amrex::FFT { /** - * \brief Poisson solver for all periodic boundaries using FFT + * \brief Poisson solver for periodic, Dirichlet & Neumann boundaries using + * FFT. */ template class Poisson @@ -40,6 +41,32 @@ private: R2X m_r2x; }; +#if (AMREX_SPACEDIM == 3) +/** + * \brief Poisson solve for Open BC using FFT. + */ +template +class PoissonOpenBC +{ +public: + + template ,int> = 0> + explicit PoissonOpenBC (Geometry const& geom, + IndexType ixtype = IndexType::TheCellType(), + IntVect const& ngrow = IntVect(0)); + + void solve (MF& soln, MF const& rhs); + + void define_doit (); // has to be public for cuda + +private: + Geometry m_geom; + Box m_grown_domain; + IntVect m_ngrow; + OpenBCSolver m_solver; +}; +#endif + /** * \brief 3D Poisson solver for periodic boundaries in the first two * dimensions and Neumann in the last dimension. @@ -123,6 +150,57 @@ void Poisson::solve (MF& soln, MF const& rhs) }); } +#if (AMREX_SPACEDIM == 3) + +template +template ,int> FOO> +PoissonOpenBC::PoissonOpenBC (Geometry const& geom, IndexType ixtype, + IntVect const& ngrow) + : m_geom(geom), + m_grown_domain(amrex::grow(amrex::convert(geom.Domain(),ixtype),ngrow)), + m_ngrow(ngrow), + m_solver(m_grown_domain) +{ + define_doit(); +} + +template +void PoissonOpenBC::define_doit () +{ + using T = typename MF::value_type; + auto const& lo = m_grown_domain.smallEnd(); + auto const dx = T(m_geom.CellSize(0)); + auto const dy = T(m_geom.CellSize(1)); + auto const dz = T(m_geom.CellSize(2)); + auto const gfac = T(1)/T(std::sqrt(T(12))); + // 0.125 comes from that there are 8 Gauss quadrature points + auto const fac = T(-0.125) * (dx*dy*dz) / (T(4)*Math::pi()); + m_solver.setGreensFunction([=] AMREX_GPU_DEVICE (int i, int j, int k) -> T + { + auto x = (T(i-lo[0]) - gfac) * dx; // first Gauss quadrature point + auto y = (T(j-lo[1]) - gfac) * dy; + auto z = (T(k-lo[2]) - gfac) * dz; + T r = 0; + for (int gx = 0; gx < 2; ++gx) { + for (int gy = 0; gy < 2; ++gy) { + for (int gz = 0; gz < 2; ++gz) { + auto xg = x + 2*gx*gfac*dx; + auto yg = y + 2*gy*gfac*dy; + auto zg = z + 2*gz*gfac*dz; + r += T(1)/std::sqrt(xg*xg+yg*yg+zg*zg); + }}} + return fac * r; + }); +} + +template +void PoissonOpenBC::solve (MF& soln, MF const& rhs) +{ + m_solver.solve(soln, rhs); +} + +#endif /* AMREX_SPACEDIM == 3 */ + template void PoissonHybrid::solve (MF& soln, MF const& rhs) { diff --git a/Src/FFT/AMReX_FFT_R2C.H b/Src/FFT/AMReX_FFT_R2C.H index 64f94b7dd81..5969c6a789b 100644 --- a/Src/FFT/AMReX_FFT_R2C.H +++ b/Src/FFT/AMReX_FFT_R2C.H @@ -10,6 +10,8 @@ namespace amrex::FFT { +template class OpenBCSolver; + /** * \brief Discrete Fourier Transform * @@ -32,6 +34,8 @@ public: MultiFab, FabArray > >; using cMF = FabArray > >; + template friend class OpenBCSolver; + /** * \brief Constructor * @@ -156,7 +160,7 @@ private: static std::pair,Plan> make_c2c_plans (cMF& inout); - void backward_doit (MF& outmf); + void backward_doit (MF& outmf, IntVect const& ngout = IntVect(0)); Plan m_fft_fwd_x{}; Plan m_fft_bwd_x{}; @@ -197,16 +201,19 @@ template R2C::R2C (Box const& domain, Info const& info) : m_real_domain(domain), m_spectral_domain_x(IntVect(0), IntVect(AMREX_D_DECL(domain.length(0)/2, - domain.bigEnd(1), - domain.bigEnd(2)))), + domain.length(1)-1, + domain.length(2)-1)), + domain.ixType()), #if (AMREX_SPACEDIM >= 2) - m_spectral_domain_y(IntVect(0), IntVect(AMREX_D_DECL(domain.bigEnd(1), + m_spectral_domain_y(IntVect(0), IntVect(AMREX_D_DECL(domain.length(1)-1, domain.length(0)/2, - domain.bigEnd(2)))), + domain.length(2)-1)), + domain.ixType()), #if (AMREX_SPACEDIM == 3) - m_spectral_domain_z(IntVect(0), IntVect(AMREX_D_DECL(domain.bigEnd(2), + m_spectral_domain_z(IntVect(0), IntVect(AMREX_D_DECL(domain.length(2)-1, domain.length(0)/2, - domain.bigEnd(1)))), + domain.length(1)-1)), + domain.ixType()), #endif #endif m_info(info) @@ -214,9 +221,7 @@ R2C::R2C (Box const& domain, Info const& info) BL_PROFILE("FFT::R2C"); static_assert(std::is_same_v || std::is_same_v); - AMREX_ALWAYS_ASSERT(m_real_domain.smallEnd() == 0 && - m_real_domain.length(0) > 1 && - m_real_domain.cellCentered()); + AMREX_ALWAYS_ASSERT(m_real_domain.length(0) > 1); #if (AMREX_SPACEDIM == 3) AMREX_ALWAYS_ASSERT(m_real_domain.length(2) > 1 || ! m_info.batch_mode); AMREX_ALWAYS_ASSERT(m_real_domain.length(1) > 1 || m_real_domain.length(2) == 1); @@ -231,13 +236,15 @@ R2C::R2C (Box const& domain, Info const& info) // make data containers // - auto bax = amrex::decompose(m_real_domain, nprocs, {AMREX_D_DECL(false,true,true)}); + auto bax = amrex::decompose(m_real_domain, nprocs, + {AMREX_D_DECL(false,true,true)}, true); DistributionMapping dmx = detail::make_iota_distromap(bax.size()); m_rx.define(bax, dmx, 1, 0, MFInfo().SetAlloc(false)); { BoxList bl = bax.boxList(); for (auto & b : bl) { + b.shift(-m_real_domain.smallEnd()); b.setBig(0, m_spectral_domain_x.bigEnd(0)); } BoxArray cbax(std::move(bl)); @@ -247,7 +254,8 @@ R2C::R2C (Box const& domain, Info const& info) #if (AMREX_SPACEDIM >= 2) DistributionMapping cdmy; if (m_real_domain.length(1) > 1) { - auto cbay = amrex::decompose(m_spectral_domain_y, nprocs, {AMREX_D_DECL(false,true,true)}); + auto cbay = amrex::decompose(m_spectral_domain_y, nprocs, + {AMREX_D_DECL(false,true,true)}, true); if (cbay.size() == dmx.size()) { cdmy = dmx; } else { @@ -261,7 +269,8 @@ R2C::R2C (Box const& domain, Info const& info) if (m_real_domain.length(1) > 1 && (! m_info.batch_mode && m_real_domain.length(2) > 1)) { - auto cbaz = amrex::decompose(m_spectral_domain_z, nprocs, {false,true,true}); + auto cbaz = amrex::decompose(m_spectral_domain_z, nprocs, + {false,true,true}, true); DistributionMapping cdmz; if (cbaz.size() == dmx.size()) { cdmz = dmx; @@ -358,7 +367,9 @@ void R2C::forward (MF const& inmf) { BL_PROFILE("FFT::R2C::forward(in)"); - m_rx.ParallelCopy(inmf, 0, 0, 1); + if (&m_rx != &inmf) { + m_rx.ParallelCopy(inmf, 0, 0, 1); + } m_fft_fwd_x.template compute_r2c(); if ( m_cmd_x2y) { @@ -380,7 +391,7 @@ void R2C::backward (MF& outmf) } template -void R2C::backward_doit (MF& outmf) +void R2C::backward_doit (MF& outmf, IntVect const& ngout) { BL_PROFILE("FFT::R2C::backward(out)"); @@ -395,7 +406,7 @@ void R2C::backward_doit (MF& outmf) } m_fft_bwd_x.template compute_r2c(); - outmf.ParallelCopy(m_rx, 0, 0, 1); + outmf.ParallelCopy(m_rx, 0, 0, 1, IntVect(0), ngout); } template diff --git a/Src/FFT/CMakeLists.txt b/Src/FFT/CMakeLists.txt index ca004337cf3..cbb205dd2e1 100644 --- a/Src/FFT/CMakeLists.txt +++ b/Src/FFT/CMakeLists.txt @@ -7,6 +7,7 @@ foreach(D IN LISTS AMReX_SPACEDIM) PRIVATE AMReX_FFT.H AMReX_FFT.cpp + AMReX_FFT_OpenBCSolver.H AMReX_FFT_R2C.H AMReX_FFT_R2X.H AMReX_FFT_Helper.H diff --git a/Src/FFT/Make.package b/Src/FFT/Make.package index 0578a8203ba..82cc4803eab 100644 --- a/Src/FFT/Make.package +++ b/Src/FFT/Make.package @@ -1,7 +1,8 @@ ifndef AMREX_FFT_MAKE AMREX_FFT_MAKE := 1 -CEXE_headers += AMReX_FFT.H AMReX_FFT_Helper.H AMReX_FFT_Poisson.H AMReX_R2C.H AMReX_R2X.H +CEXE_headers += AMReX_FFT.H AMReX_FFT_Helper.H AMReX_FFT_Poisson.H +CEXE_headers += AMReX_FFT_OpenBCSolver.H AMReX_FFT_R2C.H AMReX_FFT_R2X.H CEXE_sources += AMReX_FFT.cpp VPATH_LOCATIONS += $(AMREX_HOME)/Src/FFT diff --git a/Tests/FFT/OpenBC/CMakeLists.txt b/Tests/FFT/OpenBC/CMakeLists.txt new file mode 100644 index 00000000000..4d0690a1de5 --- /dev/null +++ b/Tests/FFT/OpenBC/CMakeLists.txt @@ -0,0 +1,13 @@ +if (NOT (3 IN_LIST AMReX_SPACEDIM) + return() +endif() + +set(_sources main.cpp) + +set(_input_files) + +setup_test(3 _sources _input_files) + +unset(_sources) +unset(_input_files) + diff --git a/Tests/FFT/OpenBC/GNUmakefile b/Tests/FFT/OpenBC/GNUmakefile new file mode 100644 index 00000000000..93376f44852 --- /dev/null +++ b/Tests/FFT/OpenBC/GNUmakefile @@ -0,0 +1,26 @@ +AMREX_HOME := ../../.. + +DEBUG = FALSE + +DIM = 3 + +COMP = gcc + +USE_MPI = TRUE +USE_OMP = FALSE +USE_CUDA = FALSE +USE_HIP = FALSE +USE_SYCL = FALSE + +USE_FFT = TRUE + +BL_NO_FORT = TRUE + +TINY_PROFILE = FALSE + +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + +include ./Make.package +include $(AMREX_HOME)/Src/Base/Make.package + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/Tests/FFT/OpenBC/Make.package b/Tests/FFT/OpenBC/Make.package new file mode 100644 index 00000000000..6b4b865e8fc --- /dev/null +++ b/Tests/FFT/OpenBC/Make.package @@ -0,0 +1 @@ +CEXE_sources += main.cpp diff --git a/Tests/FFT/OpenBC/main.cpp b/Tests/FFT/OpenBC/main.cpp new file mode 100644 index 00000000000..7049aece14a --- /dev/null +++ b/Tests/FFT/OpenBC/main.cpp @@ -0,0 +1,125 @@ +#include // Put this at the top for testing + +#include +#include +#include +#include + +#include + +using namespace amrex; + +int main (int argc, char* argv[]) +{ + static_assert(AMREX_SPACEDIM == 3); + + amrex::Initialize(argc, argv); + { + BL_PROFILE("main"); + + int n_cell_x = 128; + int n_cell_y = 128; + int n_cell_z = 128; + + int max_grid_size_x = 32; + int max_grid_size_y = 32; + int max_grid_size_z = 32; + + ParmParse pp; + pp.query("n_cell_x", n_cell_x); + pp.query("n_cell_y", n_cell_y); + pp.query("n_cell_z", n_cell_z); + pp.query("max_grid_size_x", max_grid_size_x); + pp.query("max_grid_size_y", max_grid_size_y); + pp.query("max_grid_size_z", max_grid_size_z); + + Box domain(IntVect(0), IntVect(n_cell_x-1,n_cell_y-1,n_cell_z-1)); + BoxArray ba(domain); + ba.maxSize(IntVect(max_grid_size_x, max_grid_size_y, max_grid_size_z)); + DistributionMapping dm(ba); + + Geometry geom(domain, RealBox(-1._rt, -1._rt, -1._rt, 1._rt, 1._rt, 1._rt), + CoordSys::cartesian, {AMREX_D_DECL(0,0,0)}); + + auto const& dx = geom.CellSizeArray(); + auto const& problo = geom.ProbLoArray(); + + std::array ixtypes{IndexType::TheCellType(), + IndexType::TheNodeType()}; + for (auto const ixtype : ixtypes) + { + amrex::Print() << "\nTesting " << ixtype << "\n"; + + BoxArray const& iba = amrex::convert(ba, ixtype); + int ng = ixtype.cellCentered() ? 1 : 0; + MultiFab rho(iba,dm,1,0); + MultiFab phi(iba,dm,1,ng); + phi.setVal(std::numeric_limits::max()); + + auto const& rhoma = rho.arrays(); + + constexpr int nsub = 4; + Real dxsub = dx[0]/nsub; + Real dysub = dx[1]/nsub; + Real dzsub = dx[2]/nsub; + + ParallelFor(rho, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k) + { + Real x = (i+0.5_rt/nsub)*dx[0] + problo[0]; + Real y = (j+0.5_rt/nsub)*dx[1] + problo[1]; + Real z = (k+0.5_rt/nsub)*dx[2] + problo[2]; + if (ixtype.nodeCentered()) { + x -= 0.5_rt*dx[0]; + y -= 0.5_rt*dx[1]; + z -= 0.5_rt*dx[2]; + } + int n = 0; + for (int isub = 0; isub < nsub; ++isub) { + for (int jsub = 0; jsub < nsub; ++jsub) { + for (int ksub = 0; ksub < nsub; ++ksub) { + auto xs = x + isub*dxsub; + auto ys = y + jsub*dysub; + auto zs = z + ksub*dzsub; + if ((xs*xs+ys*ys+zs*zs) < 0.25_rt) { ++n; } + }}} + rhoma[b](i,j,k) = Real(n) / Real(nsub); + }); + + FFT::PoissonOpenBC solver(geom, ixtype, IntVect(ng)); + solver.solve(phi, rho); + + Real mass = rho.sum_unique(0) * dx[0]*dx[1]*dx[2]; + Real offset = ixtype.cellCentered() ? 0.5_rt : 0.0_rt; + auto x0 = -1._rt + offset*dx[0]; + auto y0 = -1._rt + offset*dx[1]; + auto z0 = -1._rt + offset*dx[2]; + auto r0 = std::sqrt(x0*x0+y0*y0+z0*z0); // radius of the corner cell + auto expected = -mass/(4._rt*Math::pi()*r0); + amrex::Print() << " Expected phi: " << expected << "\n"; + + int iextra = ixtype.cellCentered() ? 1 : 0; + + for (int k = 0; k < 2; ++k) { + for (int j = 0; j < 2; ++j) { + for (int i = 0; i < 2; ++i) { + int ii = (i == 0) ? 0 : n_cell_x-iextra; + int jj = (j == 0) ? 0 : n_cell_y-iextra; + int kk = (k == 0) ? 0 : n_cell_z-iextra; + IntVect corner(ii,jj,kk); + auto v = amrex::get_cell_data(phi, corner); + if (!v.empty()) { + amrex::AllPrint() << " phi at " << corner << " is " << v[0] << "\n"; + auto error = std::abs(expected-v[0])/std::max(std::abs(expected),std::abs(v[0])); + amrex::AllPrint() << " error " << error << "\n"; +#ifdef AMREX_USE_FLOAT + constexpr Real eps = 1.e-5; +#else + constexpr Real eps = 1.e-6; +#endif + AMREX_ALWAYS_ASSERT(error < eps); + } + }}} + } + } + amrex::Finalize(); +}