diff --git a/Build/cmake_with_fft.sh b/Build/cmake_with_fft.sh new file mode 100755 index 000000000..0c953c7df --- /dev/null +++ b/Build/cmake_with_fft.sh @@ -0,0 +1,17 @@ + +# Example CMake config script for an OSX laptop with OpenMPI + +cmake -DCMAKE_INSTALL_PREFIX:PATH=./install \ + -DCMAKE_CXX_COMPILER:STRING=mpicxx \ + -DCMAKE_C_COMPILER:STRING=mpicc \ + -DCMAKE_Fortran_COMPILER:STRING=mpifort \ + -DMPIEXEC_PREFLAGS:STRING=--oversubscribe \ + -DCMAKE_BUILD_TYPE:STRING=Release \ + -DERF_DIM:STRING=3 \ + -DERF_ENABLE_MPI:BOOL=ON \ + -DERF_ENABLE_TESTS:BOOL=ON \ + -DERF_ENABLE_FCOMPARE:BOOL=ON \ + -DERF_ENABLE_DOCUMENTATION:BOOL=OFF \ + -DERF_ENABLE_FFT:BOOL=ON \ + -DCMAKE_EXPORT_COMPILE_COMMANDS:BOOL=ON \ + .. && make -j8 diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index f9c8a8972..884e4da91 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -47,6 +47,12 @@ function(build_erf_lib erf_lib_name) target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_EB) endif() + if(ERF_ENABLE_FFT) + target_sources(${erf_lib_name} PRIVATE + ${SRC_DIR}/Utils/ERF_solve_with_fft.cpp) + target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_FFT) + endif() + if(ERF_ENABLE_NETCDF) target_sources(${erf_lib_name} PRIVATE ${SRC_DIR}/IO/ERF_NCInterface.cpp @@ -185,6 +191,8 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Utils/ERF_MomentumToVelocity.cpp ${SRC_DIR}/Utils/ERF_PoissonSolve.cpp ${SRC_DIR}/Utils/ERF_PoissonSolve_tb.cpp + ${SRC_DIR}/Utils/ERF_solve_with_gmres.cpp + ${SRC_DIR}/Utils/ERF_solve_with_mlmg.cpp ${SRC_DIR}/Utils/ERF_TerrainMetrics.cpp ${SRC_DIR}/Utils/ERF_VelocityToMomentum.cpp ${SRC_DIR}/Utils/ERF_InteriorGhostCells.cpp diff --git a/CMake/SetAmrexOptions.cmake b/CMake/SetAmrexOptions.cmake index 99e5e78e4..038e3f0ab 100644 --- a/CMake/SetAmrexOptions.cmake +++ b/CMake/SetAmrexOptions.cmake @@ -36,6 +36,11 @@ if(ERF_ENABLE_PARTICLES) set(AMReX_PARTICLES ON) endif() +set(AMReX_FFT OFF) +if(ERF_ENABLE_FFT) + set(AMReX_FFT ON) +endif() + set(AMReX_EB OFF) if(ERF_ENABLE_EB) set(AMReX_EB ON) diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 8461f3249..b03ef317a 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -115,6 +115,10 @@ USE_LINEAR_SOLVERS_INCFLO = FALSE USE_LINEAR_SOLVERS_EM = FALSE AMReXdirs += LinearSolvers +ifeq ($(USE_FFT),TRUE) +AMReXdirs += FFT +endif + ifeq ($(USE_HDF5),TRUE) AMReXdirs += Extern/HDF5 endif diff --git a/Source/ERF.H b/Source/ERF.H index 0dcec2930..729e5d410 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -171,13 +171,20 @@ public: void project_velocities_tb (int lev, amrex::Real dt, amrex::Vector& vars, amrex::MultiFab& p); +#ifdef ERF_USE_FFT + amrex::Array,AMREX_SPACEDIM> get_fft_bc () const noexcept; + void solve_with_fft (int lev, amrex::MultiFab& rhs, amrex::MultiFab& p, + amrex::Array& fluxes); +#endif + + void solve_with_gmres (int lev, amrex::Vector& rhs, amrex::Vector& p, + amrex::Vector>& fluxes); + void solve_with_mlmg (int lev, amrex::Vector& rhs, amrex::Vector& p, + amrex::Vector>& fluxes); + // Define the projection bc's based on the domain bc types amrex::Array get_projection_bc (amrex::Orientation::Side side) const noexcept; -#ifdef ERF_USE_FFT - amrex::Array,AMREX_SPACEDIM> - get_fft_bc () const noexcept; -#endif bool projection_has_dirichlet (amrex::Array bcs) const; // Init (NOT restart or regrid) @@ -187,7 +194,8 @@ public: void restart (); // Is it time to write a plotfile or checkpoint? - bool writeNow (const amrex::Real cur_time, const amrex::Real dt, const int nstep, const int plot_int, const amrex::Real plot_per); + bool writeNow (const amrex::Real cur_time, const amrex::Real dt, const int nstep, + const int plot_int, const amrex::Real plot_per); // Called after every level 0 timestep void post_timestep (int nstep, amrex::Real time, amrex::Real dt_lev); diff --git a/Source/ERF.cpp b/Source/ERF.cpp index f1c055767..236de3b05 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1415,6 +1415,11 @@ ERF::ReadParameters () pp.query("v", verbose); pp.query("mg_v", mg_verbose); pp.query("use_fft", use_fft); +#ifndef ERF_USE_FFT + if (use_fft) { + amrex::Abort("You must build with USE_FFT in order to set use_fft = true in your inputs file"); + } +#endif // Frequency of diagnostic output pp.query("sum_interval", sum_interval); diff --git a/Source/Utils/ERF_PoissonSolve.cpp b/Source/Utils/ERF_PoissonSolve.cpp index 8a37195c9..c0f38e866 100644 --- a/Source/Utils/ERF_PoissonSolve.cpp +++ b/Source/Utils/ERF_PoissonSolve.cpp @@ -1,76 +1,8 @@ #include "ERF.H" #include "ERF_Utils.H" -#include -#include -//#include -#include -#include - using namespace amrex; -/** - * Define the domain boundary conditions for the (optional) Poisson solve - * if we want to enforce incompressibility of the initial conditions - */ - -using BCType = LinOpBCType; - -Array -ERF::get_projection_bc (Orientation::Side side) const noexcept -{ - amrex::Array r; - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - if (geom[0].isPeriodic(dir)) { - r[dir] = LinOpBCType::Periodic; - } else { - auto bc_type = domain_bc_type[Orientation(dir,side)]; - if (bc_type == "Outflow") { - r[dir] = LinOpBCType::Dirichlet; - } else - { - r[dir] = LinOpBCType::Neumann; - } - } - } - return r; -} - -#ifdef ERF_USE_FFT -Array,AMREX_SPACEDIM> -ERF::get_fft_bc () const noexcept -{ - AMREX_ALWAYS_ASSERT(geom[0].isPeriodic(0) && geom[0].isPeriodic(1)); - - Array,AMREX_SPACEDIM> r; - - for (int dir = 0; dir <= 1; dir++) { - if (geom[0].isPeriodic(dir)) { - r[dir] = std::make_pair(FFT::Boundary::periodic,FFT::Boundary::periodic); - } - } // dir - - for (OrientationIter ori; ori != nullptr; ++ori) { - const int dir = ori().coordDir(); - if (!geom[0].isPeriodic(dir) && ori().faceDir() == Orientation::low) { - auto bc_type_lo = domain_bc_type[Orientation(dir,Orientation::low)]; - auto bc_type_hi = domain_bc_type[Orientation(dir,Orientation::high)]; - if (bc_type_lo == "Outflow" && bc_type_hi == "Outflow") { - r[dir] = std::make_pair(FFT::Boundary::odd,FFT::Boundary::odd); - } else if (bc_type_lo != "Outflow" && bc_type_hi == "Outflow") { - r[dir] = std::make_pair(FFT::Boundary::even,FFT::Boundary::odd); - } else if (bc_type_lo == "Outflow" && bc_type_hi != "Outflow") { - r[dir] = std::make_pair(FFT::Boundary::odd,FFT::Boundary::even); - } else { - r[dir] = std::make_pair(FFT::Boundary::even,FFT::Boundary::even); - } - } // not periodic - } // ori - - return r; -} -#endif - /** * Project the single-level velocity field to enforce incompressibility * Note that the level may or may not be level 0. @@ -85,14 +17,6 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult auto const dom_lo = lbound(geom[lev].Domain()); auto const dom_hi = ubound(geom[lev].Domain()); - LPInfo info; - // Allow a hidden direction if the domain is one cell wide in any lateral direction - if (dom_lo.x == dom_hi.x) { - info.setHiddenDirection(0); - } else if (dom_lo.y == dom_hi.y) { - info.setHiddenDirection(1); - } - // Make sure the solver only sees the levels over which we are solving Vector ba_tmp; ba_tmp.push_back(mom_mf[Vars::cons].boxArray()); Vector dm_tmp; dm_tmp.push_back(mom_mf[Vars::cons].DistributionMap()); @@ -125,7 +49,6 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult rho0_u_const[1] = &mom_mf[IntVars::ymom]; rho0_u_const[2] = &mom_mf[IntVars::zmom]; - Real reltol = solverChoice.poisson_reltol; Real abstol = solverChoice.poisson_abstol; // **************************************************************************** @@ -169,193 +92,33 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult Print() << "Max norm of divergence before at level " << lev << " : " << rhsnorm << std::endl; } + // Initialize phi to 0 + // (It is essential that we do this in order to fill the corners; these are never + // used but the Saxpy requires the values to be initialized.) + phi[0].setVal(0.0); + if (rhsnorm <= abstol) return; Real start_step = static_cast(ParallelDescriptor::second()); -#ifdef ERF_USE_FFT // **************************************************************************** - // FFT solve + // Choose the solver and solve // **************************************************************************** - if (use_fft) - { - AMREX_ALWAYS_ASSERT(lev == 0); - // - // No terrain or stretched grids - // This calls the full 3D FFT solver with bc's set through bc_fft - // - if (!l_use_terrain) - { - if (mg_verbose > 0) { - amrex::Print() << "Using the 3D FFT solver..." << std::endl; - } - if (!m_3D_poisson) { - auto bc_fft = get_fft_bc(); - m_3D_poisson = std::make_unique>(Geom(0),bc_fft); - } - m_3D_poisson->solve(phi[lev], rhs[lev]); - - // - // Stretched grids - // This calls the hybrid 2D FFT solver + tridiagonal in z - // - // For right now we can only do this solve for periodic in the x- and y-directions - // We assume Neumann at top and bottom z-boundaries - // This will be generalized in future - // - // - } else if (l_use_terrain && SolverChoice::terrain_is_flat) - { - if (mg_verbose > 0) { amrex::Print() << "Using the hybrid FFT solver..." << std::endl; - } - if (!m_2D_poisson) { - m_2D_poisson = std::make_unique>(Geom(0)); - } - Gpu::DeviceVector stretched_dz(dom_hi.z+1, geom[lev].CellSize(2)); - m_2D_poisson->solve(phi[lev], rhs[lev], stretched_dz); - - } else { - amrex::Abort("FFT isn't appropriate for spatially varying terrain"); - } - - phi[lev].FillBoundary(geom[lev].periodicity()); - -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(phi[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Array4 const& p_arr = phi[lev].array(mfi); - - Box const& xbx = mfi.nodaltilebox(0); - const Real dx_inv = dxInv[0]; - Array4 const& fx_arr = fluxes[lev][0].array(mfi); - ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - fx_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i-1,j,k)) * dx_inv; - }); +#ifdef ERF_USE_FFT + if (use_fft) { - Box const& ybx = mfi.nodaltilebox(1); - const Real dy_inv = dxInv[1]; - Array4 const& fy_arr = fluxes[lev][1].array(mfi); - ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - fy_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i,j-1,k)) * dy_inv; - }); + solve_with_fft(lev, rhs[0], phi[0], fluxes[0]); - Box const& zbx = mfi.nodaltilebox(2); - const Real dz_inv = dxInv[2]; - Array4 const& fz_arr = fluxes[lev][2].array(mfi); - ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (k == dom_lo.z || k == dom_hi.z+1) { - fz_arr(i,j,k) = 0.0; - } else { - fz_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i,j,k-1)) * dz_inv; - } - }); - } // mfi - } else + } else #endif + if (use_gmres) { - // Initialize phi to 0 - phi[0].setVal(0.0); - - // **************************************************************************** - // GMRES solve - // **************************************************************************** - if (use_gmres) - { - amrex::Abort("GMRES isn't ready yet"); -#if 0 - MLTerrainPoisson terrpoisson(geom_tmp, ba_tmp, dm_tmp, info); - terrpoisson.setDomainBC(bclo, bchi); - terrpoisson.setMaxOrder(2); - - terrpoisson.setZPhysNd(lev, *z_phys_nd[lev]); - - if (lev > 0) { - terrpoisson.setCoarseFineBC(nullptr, ref_ratio[lev-1], LinOpBCType::Neumann); - } - terrpoisson.setLevelBC(lev, &phi[lev]); - - MLMG mlmg(terrpoisson); - GMRESMLMG gmsolver(mlmg); - gmsolver.usePrecond(false); - gmsolver.setVerbose(mg_verbose); - gmsolver.solve(phi[0], rhs[0], reltol, abstol); - - Vector phi_vec; phi_vec.resize(1); - phi_vec[0] = &phi[0]; - terrpoisson.getFluxes(GetVecOfArrOfPtrs(fluxes), phi_vec); - -#if 0 -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(phi[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Array4 const& p_arr = phi[lev].array(mfi); - - Box const& xbx = mfi.nodaltilebox(0); - const Real dx_inv = dxInv[0]; - Array4 const& fx_arr = fluxes[lev][0].array(mfi); - ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - fx_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i-1,j,k)) * dx_inv; - }); - - Box const& ybx = mfi.nodaltilebox(1); - const Real dy_inv = dxInv[1]; - Array4 const& fy_arr = fluxes[lev][1].array(mfi); - ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - fy_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i,j-1,k)) * dy_inv; - }); - - auto const dom_lo = lbound(geom[lev].Domain()); - auto const dom_hi = ubound(geom[lev].Domain()); - - Box const& zbx = mfi.nodaltilebox(2); - const Real dz_inv = dxInv[2]; - Array4 const& fz_arr = fluxes[lev][2].array(mfi); - ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (k == dom_lo.z || k == dom_hi.z+1) { - fz_arr(i,j,k) = 0.0; - } else { - fz_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i,j,k-1)) * dz_inv; - } - }); - } // mfi -#endif -#endif - - // **************************************************************************** - // Multigrid solve - // **************************************************************************** - } else { // use MLMG + solve_with_gmres(lev, rhs, phi, fluxes); - MLPoisson mlpoisson(geom_tmp, ba_tmp, dm_tmp, info); - mlpoisson.setDomainBC(bclo, bchi); - if (lev > 0) { - mlpoisson.setCoarseFineBC(nullptr, ref_ratio[lev-1], LinOpBCType::Neumann); - } - mlpoisson.setLevelBC(0, nullptr); + } else { - MLMG mlmg(mlpoisson); - int max_iter = 100; - mlmg.setMaxIter(max_iter); - - mlmg.setVerbose(mg_verbose); - mlmg.setBottomVerbose(0); - - mlmg.solve(GetVecOfPtrs(phi), - GetVecOfConstPtrs(rhs), - reltol, abstol); - mlmg.getFluxes(GetVecOfArrOfPtrs(fluxes)); - } - } // not using fft + solve_with_mlmg(lev, rhs, phi, fluxes); + } // **************************************************************************** // Subtract dt grad(phi) from the momenta (rho0u, rho0v, Omega) @@ -416,87 +179,6 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // **************************************************************************** // Update pressure variable with phi -- note that phi is dt * change in pressure // **************************************************************************** - MultiFab::Saxpy(pmf, 1.0/l_dt, phi[0],0,0,1,0); - pmf.FillBoundary(geom[lev].periodicity()); - - // **************************************************************************** - // Impose bc's on pprime - // **************************************************************************** -#ifdef _OPENMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(mom_mf[Vars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Array4 const& pp_arr = pmf.array(mfi); - Box const& bx = mfi.tilebox(); - auto const bx_lo = lbound(bx); - auto const bx_hi = ubound(bx); - if (bx_lo.x == dom_lo.x) { - if (bclo[0] == LinOpBCType::Dirichlet) { - ParallelFor(makeSlab(bx,0,dom_lo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i-1,j,k) = -pp_arr(i,j,k); - }); - } else { - ParallelFor(makeSlab(bx,0,dom_lo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i-1,j,k) = pp_arr(i,j,k); - }); - } - } - if (bx_lo.y == dom_lo.y) { - if (bclo[1] == LinOpBCType::Dirichlet) { - ParallelFor(makeSlab(bx,1,dom_lo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i,j-1,k) = -pp_arr(i,j,k); - }); - } else { - ParallelFor(makeSlab(bx,1,dom_lo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i,j-1,k) = pp_arr(i,j,k); - }); - } - } - if (bx_lo.z == dom_lo.z) { - ParallelFor(makeSlab(bx,2,dom_lo.z), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i,j,k-1) = pp_arr(i,j,k); - }); - } - if (bx_hi.x == dom_hi.x) { - if (bchi[0] == LinOpBCType::Dirichlet) { - ParallelFor(makeSlab(bx,0,dom_hi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i+1,j,k) = -pp_arr(i,j,k); - }); - } else { - ParallelFor(makeSlab(bx,0,dom_hi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i+1,j,k) = pp_arr(i,j,k); - }); - } - } - if (bx_hi.y == dom_hi.y) { - if (bchi[1] == LinOpBCType::Dirichlet) { - ParallelFor(makeSlab(bx,1,dom_hi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i,j+1,k) = -pp_arr(i,j,k); - }); - } else { - ParallelFor(makeSlab(bx,1,dom_hi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i,j+1,k) = pp_arr(i,j,k); - }); - } - } - if (bx_hi.z == dom_hi.z) { - ParallelFor(makeSlab(bx,2,dom_hi.z), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pp_arr(i,j,k+1) = pp_arr(i,j,k); - }); - } - } // mfi - - // Now overwrite with periodic fill outside domain and fine-fine fill inside + MultiFab::Saxpy(pmf, 1.0/l_dt, phi[0],0,0,1,1); pmf.FillBoundary(geom[lev].periodicity()); } diff --git a/Source/Utils/ERF_solve_with_fft.cpp b/Source/Utils/ERF_solve_with_fft.cpp new file mode 100644 index 000000000..dfb0559b3 --- /dev/null +++ b/Source/Utils/ERF_solve_with_fft.cpp @@ -0,0 +1,242 @@ +#include "ERF.H" +#include "ERF_Utils.H" + +using namespace amrex; + +#ifdef ERF_USE_FFT +Array,AMREX_SPACEDIM> +ERF::get_fft_bc () const noexcept +{ + // + // This logic only works for level 0 + // TODO: fix for level > 0 + // + Array,AMREX_SPACEDIM> r; + + for (int dir = 0; dir <= 1; dir++) { + if (geom[0].isPeriodic(dir)) { + r[dir] = std::make_pair(FFT::Boundary::periodic,FFT::Boundary::periodic); + } + } // dir + + for (OrientationIter ori; ori != nullptr; ++ori) { + const int dir = ori().coordDir(); + if (!geom[0].isPeriodic(dir) && ori().faceDir() == Orientation::low) { + auto bc_type_lo = domain_bc_type[Orientation(dir,Orientation::low)]; + auto bc_type_hi = domain_bc_type[Orientation(dir,Orientation::high)]; + if ( (bc_type_lo == "Outflow" || bc_type_lo == "Open") && + (bc_type_hi == "Outflow" || bc_type_hi == "Open") ) { + r[dir] = std::make_pair(FFT::Boundary::odd,FFT::Boundary::odd); + } else if ( (bc_type_lo != "Outflow" && bc_type_lo != "Open") && + (bc_type_hi == "Outflow" || bc_type_hi == "Outflow") ) { + r[dir] = std::make_pair(FFT::Boundary::even,FFT::Boundary::odd); + } else if ( (bc_type_lo == "Outflow" || bc_type_lo == "Open") && + (bc_type_hi != "Outflow" && bc_type_hi != "Outflow") ) { + r[dir] = std::make_pair(FFT::Boundary::odd,FFT::Boundary::even); + } else { + r[dir] = std::make_pair(FFT::Boundary::even,FFT::Boundary::even); + } + } // not periodic + } // ori + + return r; +} + +/** + * Solve the Poisson equation using FFT + * Note that the level may or may not be level 0. + */ +void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array& fluxes) +{ + BL_PROFILE("ERF::solve_with_fft()"); + + AMREX_ALWAYS_ASSERT(use_fft); + + bool l_use_terrain = SolverChoice::terrain_type != TerrainType::None; + bool use_gmres = (l_use_terrain && !SolverChoice::terrain_is_flat); + + auto const dom_lo = lbound(geom[lev].Domain()); + auto const dom_hi = ubound(geom[lev].Domain()); + + MultiFab r_hse(base_state[lev], make_alias, BaseState::r0_comp, 1); + + auto bclo = get_projection_bc(Orientation::low); + auto bchi = get_projection_bc(Orientation::high); + + // amrex::Print() << "BCLO " << bclo[0] << " " << bclo[1] << " " << bclo[2] << std::endl; + // amrex::Print() << "BCHI " << bchi[0] << " " << bchi[1] << " " << bchi[2] << std::endl; + + auto dxInv = geom[lev].InvCellSizeArray(); + + Real reltol = solverChoice.poisson_reltol; + Real abstol = solverChoice.poisson_abstol; + + // **************************************************************************** + // FFT solve + // **************************************************************************** + AMREX_ALWAYS_ASSERT(lev == 0); + // + // No terrain or stretched grids + // This calls the full 3D FFT solver with bc's set through bc_fft + // + if (!l_use_terrain) + { + if (mg_verbose > 0) { + amrex::Print() << "Using the 3D FFT solver..." << std::endl; + } + if (!m_3D_poisson) { + auto bc_fft = get_fft_bc(); + m_3D_poisson = std::make_unique>(Geom(0),bc_fft); + } + m_3D_poisson->solve(phi, rhs); + + // + // Stretched grids + // This calls the hybrid 2D FFT solver + tridiagonal in z + // + // For right now we can only do this solve for periodic in the x- and y-directions + // We assume Neumann at top and bottom z-boundaries + // This will be generalized in future + // + // + } else if (l_use_terrain && SolverChoice::terrain_is_flat) + { + if (mg_verbose > 0) { + amrex::Print() << "Using the hybrid FFT solver..." << std::endl; + } + if (!m_2D_poisson) { + m_2D_poisson = std::make_unique>(Geom(0)); + } + // m_2D_poisson->solve(phi, rhs[lev], stretched_dz_d[lev]); + + Gpu::DeviceVector stretched_dz(dom_hi.z+1, geom[lev].CellSize(2)); + m_2D_poisson->solve(phi, rhs, stretched_dz); + + } else { + amrex::Abort("FFT isn't appropriate for spatially varying terrain"); + } + + phi.FillBoundary(geom[lev].periodicity()); + + // **************************************************************************** + // Impose bc's on pprime + // **************************************************************************** +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(phi,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Array4 const& pp_arr = phi.array(mfi); + Box const& bx = mfi.tilebox(); + auto const bx_lo = lbound(bx); + auto const bx_hi = ubound(bx); + if (bx_lo.x == dom_lo.x) { + auto bc_type = domain_bc_type[Orientation(0,Orientation::low)]; + if (bc_type == "Outflow" || bc_type == "Open") { + ParallelFor(makeSlab(bx,0,dom_lo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i-1,j,k) = -pp_arr(i,j,k); + }); + } else { + ParallelFor(makeSlab(bx,0,dom_lo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i-1,j,k) = pp_arr(i,j,k); + }); + } + } + if (bx_lo.y == dom_lo.y) { + auto bc_type = domain_bc_type[Orientation(1,Orientation::low)]; + if (bc_type == "Outflow" || bc_type == "Open") { + ParallelFor(makeSlab(bx,1,dom_lo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j-1,k) = -pp_arr(i,j,k); + }); + } else { + ParallelFor(makeSlab(bx,1,dom_lo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j-1,k) = pp_arr(i,j,k); + }); + } + } + if (bx_lo.z == dom_lo.z) { + ParallelFor(makeSlab(bx,2,dom_lo.z), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j,k-1) = pp_arr(i,j,k); + }); + } + if (bx_hi.x == dom_hi.x) { + auto bc_type = domain_bc_type[Orientation(0,Orientation::high)]; + if (bc_type == "Outflow" || bc_type == "Open") { + ParallelFor(makeSlab(bx,0,dom_hi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i+1,j,k) = -pp_arr(i,j,k); + }); + } else { + ParallelFor(makeSlab(bx,0,dom_hi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i+1,j,k) = pp_arr(i,j,k); + }); + } + } + if (bx_hi.y == dom_hi.y) { + auto bc_type = domain_bc_type[Orientation(1,Orientation::high)]; + if (bc_type == "Outflow" || bc_type == "Open") { + ParallelFor(makeSlab(bx,1,dom_hi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j+1,k) = -pp_arr(i,j,k); + }); + } else { + ParallelFor(makeSlab(bx,1,dom_hi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j+1,k) = pp_arr(i,j,k); + }); + } + } + if (bx_hi.z == dom_hi.z) { + ParallelFor(makeSlab(bx,2,dom_hi.z), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j,k+1) = pp_arr(i,j,k); + }); + } + } // mfi + + // Now overwrite with periodic fill outside domain and fine-fine fill inside + phi.FillBoundary(geom[lev].periodicity()); + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(phi, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Array4 const& p_arr = phi.array(mfi); + + Box const& xbx = mfi.nodaltilebox(0); + const Real dx_inv = dxInv[0]; + Array4 const& fx_arr = fluxes[0].array(mfi); + ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + fx_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i-1,j,k)) * dx_inv; + }); + + Box const& ybx = mfi.nodaltilebox(1); + const Real dy_inv = dxInv[1]; + Array4 const& fy_arr = fluxes[1].array(mfi); + ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + fy_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i,j-1,k)) * dy_inv; + }); + + Box const& zbx = mfi.nodaltilebox(2); + const Real dz_inv = dxInv[2]; + Array4 const& fz_arr = fluxes[2].array(mfi); + ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (k == dom_lo.z || k == dom_hi.z+1) { + fz_arr(i,j,k) = 0.0; + } else { + fz_arr(i,j,k) = -(p_arr(i,j,k) - p_arr(i,j,k-1)) * dz_inv; + } + }); + } // mfi +} +#endif diff --git a/Source/Utils/ERF_solve_with_gmres.cpp b/Source/Utils/ERF_solve_with_gmres.cpp new file mode 100644 index 000000000..cd2d3465c --- /dev/null +++ b/Source/Utils/ERF_solve_with_gmres.cpp @@ -0,0 +1,66 @@ +#include "ERF.H" +#include "ERF_Utils.H" + +#include +#include +#include +#include + +using namespace amrex; + +/** + * Solve the Poisson equation using MLMG + * Note that the level may or may not be level 0. + */ +void ERF::solve_with_gmres (int lev, Vector& /*rhs*/, Vector& /*phi*/, Vector>& /*fluxes*/) +{ + BL_PROFILE("ERF::solve_with_mlmg()"); + + auto const dom_lo = lbound(geom[lev].Domain()); + auto const dom_hi = ubound(geom[lev].Domain()); + + LPInfo info; + // Allow a hidden direction if the domain is one cell wide in any lateral direction + if (dom_lo.x == dom_hi.x) { + info.setHiddenDirection(0); + } else if (dom_lo.y == dom_hi.y) { + info.setHiddenDirection(1); + } + + // Make sure the solver only sees the levels over which we are solving + Vector ba_tmp; ba_tmp.push_back(rhs[lev].boxArray()); + Vector dm_tmp; dm_tmp.push_back(rhs[lev].DistributionMap()); + Vector geom_tmp; geom_tmp.push_back(geom[lev]); + + auto bclo = get_projection_bc(Orientation::low); + auto bchi = get_projection_bc(Orientation::high); + + // amrex::Print() << "BCLO " << bclo[0] << " " << bclo[1] << " " << bclo[2] << std::endl; + // amrex::Print() << "BCHI " << bchi[0] << " " << bchi[1] << " " << bchi[2] << std::endl; + +#if 0 + Real reltol = solverChoice.poisson_reltol; + Real abstol = solverChoice.poisson_abstol; + + MLTerrainPoisson terrpoisson(geom_tmp, ba_tmp, dm_tmp, info); + terrpoisson.setDomainBC(bclo, bchi); + terrpoisson.setMaxOrder(2); + + terrpoisson.setZPhysNd(lev, *z_phys_nd[lev]); + + if (lev > 0) { + terrpoisson.setCoarseFineBC(nullptr, ref_ratio[lev-1], LinOpBCType::Neumann); + } + terrpoisson.setLevelBC(lev, &phi[lev]); + + MLMG mlmg(terrpoisson); + GMRESMLMG gmsolver(mlmg); + gmsolver.usePrecond(false); + gmsolver.setVerbose(mg_verbose); + gmsolver.solve(phi[0], rhs[0], reltol, abstol); + + Vector phi_vec; phi_vec.resize(1); + phi_vec[0] = &phi[0]; + terrpoisson.getFluxes(GetVecOfArrOfPtrs(fluxes), phi_vec); +#endif +} diff --git a/Source/Utils/ERF_solve_with_mlmg.cpp b/Source/Utils/ERF_solve_with_mlmg.cpp new file mode 100644 index 000000000..49e4ef532 --- /dev/null +++ b/Source/Utils/ERF_solve_with_mlmg.cpp @@ -0,0 +1,91 @@ +#include "ERF.H" +#include "ERF_Utils.H" + +#include +#include + +using namespace amrex; + +/** + * Define the domain boundary conditions for the (optional) Poisson solve + * if we want to enforce incompressibility of the initial conditions + */ + +using BCType = LinOpBCType; + +Array +ERF::get_projection_bc (Orientation::Side side) const noexcept +{ + amrex::Array r; + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + if (geom[0].isPeriodic(dir)) { + r[dir] = LinOpBCType::Periodic; + } else { + auto bc_type = domain_bc_type[Orientation(dir,side)]; + if (bc_type == "Outflow") { + r[dir] = LinOpBCType::Dirichlet; + } else + { + r[dir] = LinOpBCType::Neumann; + } + } + } + return r; +} + +/** + * Solve the Poisson equation using MLMG + * Note that the level may or may not be level 0. + */ +void ERF::solve_with_mlmg (int lev, Vector& rhs, Vector& phi, Vector>& fluxes) +{ + BL_PROFILE("ERF::solve_with_mlmg()"); + + auto const dom_lo = lbound(geom[lev].Domain()); + auto const dom_hi = ubound(geom[lev].Domain()); + + LPInfo info; + // Allow a hidden direction if the domain is one cell wide in any lateral direction + if (dom_lo.x == dom_hi.x) { + info.setHiddenDirection(0); + } else if (dom_lo.y == dom_hi.y) { + info.setHiddenDirection(1); + } + + // Make sure the solver only sees the levels over which we are solving + Vector ba_tmp; ba_tmp.push_back(rhs[lev].boxArray()); + Vector dm_tmp; dm_tmp.push_back(rhs[lev].DistributionMap()); + Vector geom_tmp; geom_tmp.push_back(geom[lev]); + + auto bclo = get_projection_bc(Orientation::low); + auto bchi = get_projection_bc(Orientation::high); + + // amrex::Print() << "BCLO " << bclo[0] << " " << bclo[1] << " " << bclo[2] << std::endl; + // amrex::Print() << "BCHI " << bchi[0] << " " << bchi[1] << " " << bchi[2] << std::endl; + + Real reltol = solverChoice.poisson_reltol; + Real abstol = solverChoice.poisson_abstol; + + // **************************************************************************** + // Multigrid solve + // **************************************************************************** + + MLPoisson mlpoisson(geom_tmp, ba_tmp, dm_tmp, info); + mlpoisson.setDomainBC(bclo, bchi); + if (lev > 0) { + mlpoisson.setCoarseFineBC(nullptr, ref_ratio[lev-1], LinOpBCType::Neumann); + } + mlpoisson.setLevelBC(0, nullptr); + + MLMG mlmg(mlpoisson); + int max_iter = 100; + mlmg.setMaxIter(max_iter); + + mlmg.setVerbose(mg_verbose); + mlmg.setBottomVerbose(0); + + mlmg.solve(GetVecOfPtrs(phi), + GetVecOfConstPtrs(rhs), + reltol, abstol); + mlmg.getFluxes(GetVecOfArrOfPtrs(fluxes)); +} diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index 238eccb8c..3fc5cd0ec 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -26,3 +26,6 @@ CEXE_sources += ERF_Time_Avg_Vel.cpp CEXE_sources += ERF_PoissonSolve.cpp CEXE_sources += ERF_PoissonSolve_tb.cpp +CEXE_sources += ERF_solve_with_fft.cpp +CEXE_sources += ERF_solve_with_gmres.cpp +CEXE_sources += ERF_solve_with_mlmg.cpp diff --git a/Submodules/AMReX b/Submodules/AMReX index 4a6e7c87e..456c93c7d 160000 --- a/Submodules/AMReX +++ b/Submodules/AMReX @@ -1 +1 @@ -Subproject commit 4a6e7c87e1494d10893b36f0d7ce6bf9df19fe0d +Subproject commit 456c93c7d9512f1cdffac0574973d7df41417898