diff --git a/Projections/CMakeLists.txt b/Projections/CMakeLists.txt index a2c4c5d0f..773ebce18 100644 --- a/Projections/CMakeLists.txt +++ b/Projections/CMakeLists.txt @@ -12,3 +12,9 @@ target_sources( hydro_NodalProjector.cpp hydro_NodalProjector.H ) + +if (HYDRO_FFT) + target_sources(amrex_hydro PRIVATE + hydro_FFTMacProjector.cpp + hydro_FFTMacProjector.H) +endif() diff --git a/Projections/Make.package b/Projections/Make.package index 339610577..b398b3870 100644 --- a/Projections/Make.package +++ b/Projections/Make.package @@ -3,3 +3,8 @@ CEXE_headers += hydro_NodalProjector.H CEXE_sources += hydro_MacProjector.cpp CEXE_sources += hydro_NodalProjector.cpp + +ifeq ($(USE_FFT),TRUE) + CEXE_sources += hydro_FFTMacProjector.cpp + CEXE_headers += hydro_FFTMacProjector.H +endif diff --git a/Projections/hydro_FFTMacProjector.H b/Projections/hydro_FFTMacProjector.H new file mode 100644 index 000000000..3ec7e1fc7 --- /dev/null +++ b/Projections/hydro_FFTMacProjector.H @@ -0,0 +1,30 @@ +#ifndef AMREX_FFT_MAC_PROJECTOR_H_ +#define AMREX_FFT_MAC_PROJECTOR_H_ + +#include +#include + +namespace Hydro { + +class FFTMacProjector +{ +public: + + FFTMacProjector (amrex::Geometry const& geom, + amrex::Array const& lobc, + amrex::Array const& hibc); + + void setUMAC (amrex::Array const& umac); + + void project (); + +private: + amrex::Geometry m_geom; + amrex::Array m_lobc; + amrex::Array m_hibc; + amrex::FFT::Poisson m_fft; + amrex::Array m_umac; +}; + +} +#endif diff --git a/Projections/hydro_FFTMacProjector.cpp b/Projections/hydro_FFTMacProjector.cpp new file mode 100644 index 000000000..253edbb2f --- /dev/null +++ b/Projections/hydro_FFTMacProjector.cpp @@ -0,0 +1,114 @@ +#include +#include + +using namespace amrex; + +namespace Hydro { + +namespace { + Array,AMREX_SPACEDIM> + make_fft_bc (amrex::Array const& lobc, + amrex::Array const& hibc) + { + Array,AMREX_SPACEDIM> fft_bc; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if ((lobc[idim] == LinOpBCType::Periodic && + hibc[idim] != LinOpBCType::Periodic) || + (lobc[idim] != LinOpBCType::Periodic && + hibc[idim] == LinOpBCType::Periodic)) + { + amrex::Abort("FFTMacProjector: wrong periodic BC type"); + } + for (int face = 0; face < 2; ++face) { + auto& lbc = (face == 0) ? lobc[idim] : hibc[idim]; + auto& fbc = (face == 0) ? fft_bc[idim].first : fft_bc[idim].second; + switch (lbc) + { + case LinOpBCType::Periodic: + fbc = FFT::Boundary::periodic; + break; + case LinOpBCType::Dirichlet: + fbc = FFT::Boundary::odd; + break; + case LinOpBCType::Neumann: + fbc = FFT::Boundary::even; + break; + default: + amrex::Abort("FFTMacProjector: only Periodic, Dirichlet & Neumann are supported"); + } + } + } + return fft_bc; + } +} + +FFTMacProjector::FFTMacProjector ( + amrex::Geometry const& geom, + amrex::Array const& lobc, + amrex::Array const& hibc) + : m_geom(geom), + m_lobc(lobc), + m_hibc(hibc), + m_fft(geom, make_fft_bc(lobc,hibc)) +{} + +void FFTMacProjector::setUMAC (Array const& umac) +{ + m_umac = umac; +} + +void FFTMacProjector::project () +{ + MultiFab rhs(amrex::convert(m_umac[0]->boxArray(), IntVect(0)), + m_umac[0]->DistributionMap(), 1, 0); + amrex::computeDivergence(rhs, GetArrOfConstPtrs(m_umac), m_geom); + + AMREX_ALWAYS_ASSERT(m_geom.Domain().numPts() == rhs.boxArray().numPts()); + + bool has_dirichlet = false; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + has_dirichlet = has_dirichlet || + m_lobc[idim] == LinOpBCType::Dirichlet || + m_hibc[idim] == LinOpBCType::Dirichlet; + } + if (! has_dirichlet) { + auto rhosum = rhs.sum(0); + rhs.plus(-rhosum/m_geom.Domain().d_numPts(), 0, 1); + } + + MultiFab phi(rhs.boxArray(), rhs.DistributionMap(), 1, 1); + m_fft.solve(phi, rhs); + phi.FillBoundary(m_geom.periodicity()); + + auto const& phima = phi.const_arrays(); + Box const& domain = m_geom.growPeriodicDomain(1); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + auto const& uma = m_umac[idim]->arrays(); + int domlo = domain.smallEnd(idim); + int domhi = domain.bigEnd(idim)+1; + auto const bclo = m_lobc[idim]; + auto const bchi = m_hibc[idim]; + Real const dxinv = m_geom.InvCellSize(idim); + ParallelFor(*m_umac[idim], [=] AMREX_GPU_DEVICE (int b, int i, int j, int k) + { + IntVect iv(AMREX_D_DECL(i,j,k)); + IntVect miv = iv - IntVect::TheDimensionVector(idim); + auto const& u = uma[b]; + auto const& p = phima[b]; + if (iv[idim] == domlo) { + if (bclo == LinOpBCType::Dirichlet) { + u(i,j,k) -= Real(2.0)*dxinv*p(iv); + } + } else if (iv[idim] == domhi) { + if (bchi == LinOpBCType::Dirichlet) { + u(i,j,k) += Real(2.0)*dxinv*p(miv); + } + } else { + u(i,j,k) -= dxinv*(p(iv)-p(miv)); + } + }); + } + Gpu::streamSynchronize(); +} + +} diff --git a/Projections/hydro_MacProjector.H b/Projections/hydro_MacProjector.H index d9c789e7d..1272f4a3b 100644 --- a/Projections/hydro_MacProjector.H +++ b/Projections/hydro_MacProjector.H @@ -10,14 +10,14 @@ #include #endif -#ifdef AMREX_USE_FFT -#include -#endif - #ifdef AMREX_USE_HYPRE #include #endif +#ifdef AMREX_USE_FFT +#include +#endif + namespace Hydro { class MacProjector @@ -216,12 +216,6 @@ private: #ifdef AMREX_USE_HYPRE std::unique_ptr m_hypremlabeclap; #endif - - bool m_use_fft = false; -#ifdef AMREX_USE_FFT - std::unique_ptr> m_fft_poisson; - std::unique_ptr> m_fft_poisson_hybrid; -#endif }; } diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index 391d32207..d7da86c69 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -24,7 +24,7 @@ MacProjector::MacProjector( m_phi_loc(a_phi_loc), m_divu_loc(a_divu_loc) { - amrex::ignore_unused(m_divu_loc, m_beta_loc, m_phi_loc, m_umac_loc, m_use_fft); + amrex::ignore_unused(m_divu_loc, m_beta_loc, m_phi_loc, m_umac_loc); } MacProjector::MacProjector (const Vector >& a_umac, @@ -259,13 +259,6 @@ MacProjector::setDomainBC (const Array& lobc, #ifdef AMREX_USE_HYPRE m_hypremlabeclap.reset(); #endif -#ifdef AMREX_USE_FFT - if (m_fft_poisson_hybrid) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lobc.back() == LinOpBCType::Neumann && - hibc.back() == LinOpBCType::Neumann, - "FFT::PoissonHybrid supports Neumann BC in z-direction only"); - } -#endif } @@ -394,18 +387,6 @@ MacProjector::project (Real reltol, Real atol) m_mlmg->prepareForFluxes(amrex::GetVecOfConstPtrs(m_phi)); } } else -#endif -#ifdef AMREX_USE_FFT - if (m_fft_poisson || m_fft_poisson_hybrid) { - m_mlmg->makeSolvable(0,0,m_rhs[0]); - if (m_fft_poisson) { - m_fft_poisson->solve(m_phi[0], m_rhs[0]); - } else { - m_fft_poisson_hybrid->solve(m_phi[0], m_rhs[0]); - } - // The solve below should not do any iterations. - m_mlmg->solve(amrex::GetVecOfPtrs(m_phi), amrex::GetVecOfConstPtrs(m_rhs), reltol, atol); - } else #endif { m_mlmg->solve(amrex::GetVecOfPtrs(m_phi), amrex::GetVecOfConstPtrs(m_rhs), reltol, atol); @@ -584,27 +565,6 @@ void MacProjector::initProjector (Vector const& a_grids, } auto const& dm = a_dmap; -#ifdef AMREX_USE_FFT - { - ParmParse pp("mac_proj"); - pp.query("use_fft", m_use_fft); - } - if (m_use_fft) { - bool fft_supported = (nlevs == 1) && (a_overset_mask.empty()) && - m_geom[0].Domain().numPts() == ba[0].numPts(); -#if (AMREX_SPACEDIM == 3) - fft_supported = fft_supported && (m_geom[0].isPeriodic(0) && - m_geom[0].isPeriodic(1)); -#else - fft_supported = fft_supported && m_geom[0].isAllPeriodic(); -#endif - if (!fft_supported) { - m_use_fft = false; - amrex::Warning("MacProjector: FFT not support"); - } - } -#endif - m_rhs.resize(nlevs); m_phi.resize(nlevs); m_fluxes.resize(nlevs); @@ -622,7 +582,7 @@ void MacProjector::initProjector (Vector const& a_grids, } } - if (m_use_mlhypre || m_use_fft) { a_lpinfo.setMaxCoarseningLevel(0); } + if (m_use_mlhypre) { a_lpinfo.setMaxCoarseningLevel(0); } if (a_overset_mask.empty()) { m_poisson = std::make_unique(m_geom, ba, dm, a_lpinfo); @@ -648,16 +608,6 @@ void MacProjector::initProjector (Vector const& a_grids, } #endif -#ifdef AMREX_USE_FFT - if (m_use_fft) { - if (m_geom[0].isAllPeriodic()) { - m_fft_poisson = std::make_unique>(m_geom[0]); - } else { - m_fft_poisson_hybrid = std::make_unique>(m_geom[0]); - } - } -#endif - setOptions(); m_needs_init = false; diff --git a/Tests/FFT_Mac_Projection/GNUmakefile b/Tests/FFT_Mac_Projection/GNUmakefile new file mode 100644 index 000000000..5d155a18f --- /dev/null +++ b/Tests/FFT_Mac_Projection/GNUmakefile @@ -0,0 +1,41 @@ +AMREX_HOME ?= ../../../amrex +AMREX_HYDRO_HOME = ../.. + +USE_FFT = TRUE +USE_MPI = TRUE +USE_OMP = FALSE + +COMP = gnu + +DIM = 3 + +DEBUG = FALSE + +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + +include ./Make.package + +Pdirs := AmrCore +Pdirs += Base +Pdirs += Boundary +Pdirs += FFT +Pdirs += LinearSolvers/MLMG + +Ppack += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package) +Ppack += $(AMREX_HYDRO_HOME)/Projections/Make.package + +include $(Ppack) + +Bdirs := AmrCore +Bdirs += Base +Bdirs += Boundary +Bdirs += FFT +Bdirs += EB + +Blocs := $(foreach dir, $(Bdirs), $(AMREX_HOME)/Src/$(dir)) +Blocs += $(AMREX_HYDRO_HOME)/Projections + +INCLUDE_LOCATIONS += $(Blocs) +VPATH_LOCATIONS += $(Blocs) + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/Tests/FFT_Mac_Projection/Make.package b/Tests/FFT_Mac_Projection/Make.package new file mode 100644 index 000000000..6b4b865e8 --- /dev/null +++ b/Tests/FFT_Mac_Projection/Make.package @@ -0,0 +1 @@ +CEXE_sources += main.cpp diff --git a/Tests/FFT_Mac_Projection/main.cpp b/Tests/FFT_Mac_Projection/main.cpp new file mode 100644 index 000000000..d10c605c6 --- /dev/null +++ b/Tests/FFT_Mac_Projection/main.cpp @@ -0,0 +1,71 @@ +#include +#include +#include + +#include + +using namespace amrex; +using namespace Hydro; + +int main (int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + { + int n_cell = 128; + int max_grid_size = 32; + { + ParmParse pp; + pp.query("n_cell", n_cell); + pp.query("max_grid_size", max_grid_size); + } + + Geometry geom; + BoxArray grids; + DistributionMapping dmap; + { + RealBox rb({AMREX_D_DECL(0.,0.,0.)}, {AMREX_D_DECL(1.,1.,1.)}); + + Array isp{AMREX_D_DECL(0,1,0)}; + Box domain(IntVect(0), IntVect(n_cell-1)); + geom.define(domain, rb, 0, isp); + + grids.define(domain); + grids.maxSize(max_grid_size); + + dmap.define(grids); + } + + Array vel; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + vel[idim].define(amrex::convert(grids,IntVect::TheDimensionVector(idim)), + dmap, 1, 0); + vel[idim].setVal(0); + Box region = geom.Domain(); + region.growLo(idim, -geom.Domain().length(idim)/4); + region.growHi(idim, -geom.Domain().length(idim)/4); + vel[idim].setVal(Real(1.0), region, 0, 1); + } + + MultiFab divu(grids, dmap, 1, 0); + amrex::computeDivergence(divu, GetArrOfConstPtrs(vel), geom); + amrex::Print() << "Divergence before projection: " + << divu.min(0) << " " << divu.max(0) << "\n"; + + FFTMacProjector macproj(geom, + {AMREX_D_DECL(LinOpBCType::Dirichlet, + LinOpBCType::Periodic, + LinOpBCType::Neumann)}, + {AMREX_D_DECL(LinOpBCType::Neumann, + LinOpBCType::Periodic, + LinOpBCType::Dirichlet)}); + + macproj.setUMAC(GetArrOfPtrs(vel)); + + macproj.project(); + + amrex::computeDivergence(divu, GetArrOfConstPtrs(vel), geom); + amrex::Print() << "Divergence after projection: " + << divu.min(0) << " " << divu.max(0) << "\n"; + } + amrex::Finalize(); +}