diff --git a/Src/Base/AMReX_FabDataType.H b/Src/Base/AMReX_FabDataType.H new file mode 100644 index 00000000000..81537ae8065 --- /dev/null +++ b/Src/Base/AMReX_FabDataType.H @@ -0,0 +1,27 @@ +#ifndef AMREX_FAB_DATA_TYPE_H_ +#define AMREX_FAB_DATA_TYPE_H_ +#include + +#include + +namespace amrex { + +template struct FabDataType {}; +// +template +struct FabDataType > > +{ + using fab_type = typename T::fab_type; + using value_type = typename T::value_type; +}; + +template +struct FabDataType > > +{ + using fab_type = typename T::value_type::fab_type; + using value_type = typename T::value_type::value_type; +}; + +} + +#endif diff --git a/Src/Base/CMakeLists.txt b/Src/Base/CMakeLists.txt index 459ec3bd7c4..ab30abc10f0 100644 --- a/Src/Base/CMakeLists.txt +++ b/Src/Base/CMakeLists.txt @@ -132,6 +132,7 @@ foreach(D IN LISTS AMReX_SPACEDIM) AMReX_Array4.H AMReX_MakeType.H AMReX_TypeTraits.H + AMReX_FabDataType.H AMReX_FabFactory.H AMReX_BaseFabUtility.H # Fortran data defined on unions of rectangles ---------------------------- diff --git a/Src/Base/Make.package b/Src/Base/Make.package index 276887ebd79..ed0293a957a 100644 --- a/Src/Base/Make.package +++ b/Src/Base/Make.package @@ -149,6 +149,7 @@ C$(AMREX_BASE)_headers += AMReX_IArrayBox.H C$(AMREX_BASE)_headers += AMReX_MakeType.H C$(AMREX_BASE)_headers += AMReX_TypeTraits.H +C$(AMREX_BASE)_headers += AMReX_FabDataType.H C$(AMREX_BASE)_headers += AMReX_Array4.H C$(AMREX_BASE)_sources += AMReX_BaseFab.cpp diff --git a/Src/Boundary/AMReX_FabSet.H b/Src/Boundary/AMReX_FabSet.H index 9841555b336..0e9051b182d 100644 --- a/Src/Boundary/AMReX_FabSet.H +++ b/Src/Boundary/AMReX_FabSet.H @@ -3,6 +3,7 @@ #define AMREX_FABSET_H_ #include +#include #include #include #include @@ -46,8 +47,8 @@ class FabSetT friend class FabSetIter; friend class FluxRegister; public: - using value_type = typename MF::value_type; - using FAB = typename MF::fab_type; + using value_type = typename FabDataType::value_type; + using FAB = typename FabDataType::fab_type; // //! The default constructor -- you must later call define(). diff --git a/Src/Boundary/AMReX_LO_BCTYPES.H b/Src/Boundary/AMReX_LO_BCTYPES.H index 2b4e9178bb2..e8fdb3a0274 100644 --- a/Src/Boundary/AMReX_LO_BCTYPES.H +++ b/Src/Boundary/AMReX_LO_BCTYPES.H @@ -20,6 +20,7 @@ #define AMREX_LO_INFLOW 106 #define AMREX_LO_INHOMOG_NEUMANN 107 #define AMREX_LO_ROBIN 108 +#define AMREX_LO_SYMMETRY 109 #define AMREX_LO_PERIODIC 200 #define AMREX_LO_BOGUS 1729 @@ -38,6 +39,7 @@ namespace amrex { inflow = AMREX_LO_INFLOW, inhomogNeumann = AMREX_LO_INHOMOG_NEUMANN, Robin = AMREX_LO_ROBIN, + symmetry = AMREX_LO_SYMMETRY, Periodic = AMREX_LO_PERIODIC, bogus = AMREX_LO_BOGUS }; @@ -48,4 +50,3 @@ namespace amrex { #endif #endif - diff --git a/Src/Boundary/AMReX_LO_BCTYPES.cpp b/Src/Boundary/AMReX_LO_BCTYPES.cpp index 85731b6af47..f9f78c766bb 100644 --- a/Src/Boundary/AMReX_LO_BCTYPES.cpp +++ b/Src/Boundary/AMReX_LO_BCTYPES.cpp @@ -52,6 +52,11 @@ std::ostream& operator<< (std::ostream& os, const LinOpBCType& t) os << "Robin"; break; } + case LinOpBCType::symmetry: + { + os << "symmetry"; + break; + } case LinOpBCType::Periodic: { os << "Periodic"; diff --git a/Src/LinearSolvers/CMakeLists.txt b/Src/LinearSolvers/CMakeLists.txt index 76f75e06123..cae0b2028f0 100644 --- a/Src/LinearSolvers/CMakeLists.txt +++ b/Src/LinearSolvers/CMakeLists.txt @@ -67,6 +67,15 @@ foreach(D IN LISTS AMReX_SPACEDIM) ) endif () + if (NOT D EQUAL 1) + target_sources(amrex_${D}d + PRIVATE + MLMG/AMReX_MLCurlCurl.H + MLMG/AMReX_MLCurlCurl.cpp + MLMG/AMReX_MLCurlCurl_K.H + ) + endif () + if (AMReX_EB AND NOT D EQUAL 1) target_sources(amrex_${D}d PRIVATE diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H new file mode 100644 index 00000000000..a366bb65567 --- /dev/null +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H @@ -0,0 +1,155 @@ +#ifndef AMREX_ML_CURL_CURL_H_ +#define AMREX_ML_CURL_CURL_H_ +#include + +#include + +namespace amrex { + +/** + * \brief curl (alpha curl E) + beta E = rhs + * + * Here E is an Array of 3 MultiFabs on staggered grid, alpha is a positive + * scalar, and beta is a non-negative scalar. + * + * TODO: If beta is zero, the system could be singular. + * TODO: Try different restriction & interpolation strategies. + */ +class MLCurlCurl + : public MLLinOpT > +{ +public: + using MF = Array; + using RT = typename MLLinOpT::RT; + using BCType = typename MLLinOpT::BCType; + using BCMode = typename MLLinOpT::BCMode; + using StateMode = typename MLLinOpT::StateMode; + using Location = typename MLLinOpT::Location; + + MLCurlCurl () = default; + MLCurlCurl (const Vector& a_geom, + const Vector& a_grids, + const Vector& a_dmap, + const LPInfo& a_info = LPInfo()); + + void define (const Vector& a_geom, + const Vector& a_grids, + const Vector& a_dmap, + const LPInfo& a_info = LPInfo()); + + void setScalars (RT a_alpha, RT a_beta) noexcept; + + [[nodiscard]] std::string name () const override { + return std::string("curl of curl"); + } + + void setLevelBC (int amrlev, const MF* levelbcdata, + const MF* robinbc_a = nullptr, + const MF* robinbc_b = nullptr, + const MF* robinbc_f = nullptr) override; + + void restriction (int amrlev, int cmglev, MF& crse, MF& fine) const override; + + void interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const override; + + void interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const override; + + void interpolationAmr (int famrlev, MF& fine, const MF& crse, + IntVect const& nghost) const override; + + void averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs, + const MF& fine_sol, const MF& fine_rhs) override; + + void apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode, + StateMode s_mode, const MLMGBndryT* bndry=nullptr) const override; + + void smooth (int amrlev, int mglev, MF& sol, const MF& rhs, + bool skip_fillboundary=false) const override; + + void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, + const MF* crse_bcdata=nullptr) override; + + void correctionResidual (int amrlev, int mglev, MF& resid, MF& x, + const MF& b, BCMode bc_mode, + const MF* crse_bcdata=nullptr) override; + + void reflux (int crse_amrlev, + MF& res, const MF& crse_sol, const MF& crse_rhs, + MF& fine_res, MF& fine_sol, const MF& fine_rhs) const override; + + void compFlux (int amrlev, const Array& fluxes, + MF& sol, Location loc) const override; + + void compGrad (int amrlev, const Array& grad, + MF& sol, Location loc) const override; + + void applyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const override {} + + void unapplyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const override {} + + void prepareForSolve () override; + + [[nodiscard]] bool isSingular (int /*amrlev*/) const override { return false; } + [[nodiscard]] bool isBottomSingular () const override { return false; } + + RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const override; + + std::unique_ptr> makeNLinOp (int /*grid_size*/) const override { return nullptr; } + + [[nodiscard]] RT normInf (int amrlev, MF const& mf, bool local) const override; + + void averageDownAndSync (Vector& sol) const override; + + void avgDownResAmr (int clev, MF& cres, MF const& fres) const override; + + void avgDownResMG (int clev, MF& cres, MF const& fres) const override; + + [[nodiscard]] IntVect getNGrowVectRestriction () const override { + return IntVect(0); + } + + void make (Vector >& mf, IntVect const& ng) const override; + + [[nodiscard]] MF make (int amrlev, int mglev, IntVect const& ng) const override; + + [[nodiscard]] MF makeAlias (MF const& mf) const override; + + [[nodiscard]] MF makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const override; + + [[nodiscard]] MF makeCoarseAmr (int famrlev, IntVect const& ng) const override; + +// public for cuda + + void smooth (int amrlev, int mglev, MF& sol, MultiFab const& rhs, + int redblack) const; + + void compresid (int amrlev, int mglev, MF& resid, MF const& b) const; + + void applyPhysBC (int amrlev, int mglev, MultiFab& mf) const; + +private: + + void applyBC (int amrlev, int mglev, MF& in) const; + void applyBC (int amrlev, int mglev, MultiFab& mf) const; + + [[nodiscard]] iMultiFab const& getDotMask (int amrlev, int mglev, int idim) const; + + [[nodiscard]] int getDirichlet (int amrlev, int mglev, int idim, int face) const; + + RT m_alpha = std::numeric_limits::lowest(); + RT m_beta = std::numeric_limits::lowest(); + + Array m_etype +#if (AMREX_SPACEDIM == 3) + {IntVect(0,1,1), IntVect(1,0,1), IntVect(1,1,0)}; +#else + {IntVect(0,1), IntVect(1,0), IntVect(1,1)}; +#endif + + mutable Vector,3>>> m_dotmask; + static constexpr int m_ncomp = 1; +}; + +} + +#endif diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp new file mode 100644 index 00000000000..94d2b1d17ce --- /dev/null +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp @@ -0,0 +1,552 @@ + +#include +#include + +namespace amrex { + +MLCurlCurl::MLCurlCurl (const Vector& a_geom, + const Vector& a_grids, + const Vector& a_dmap, + const LPInfo& a_info) +{ + define(a_geom, a_grids, a_dmap, a_info); +} + +void MLCurlCurl::define (const Vector& a_geom, + const Vector& a_grids, + const Vector& a_dmap, + const LPInfo& a_info) +{ + MLLinOpT::define(a_geom, a_grids, a_dmap, a_info, {}); + + m_dotmask.resize(this->m_num_amr_levels); + for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) { + m_dotmask[amrlev].resize(this->m_num_mg_levels[amrlev]); + } +} + +void MLCurlCurl::setScalars (RT a_alpha, RT a_beta) noexcept +{ + m_alpha = a_alpha; + m_beta = a_beta; +} + +void MLCurlCurl::setLevelBC (int amrlev, const MF* levelbcdata, // TODO + const MF* robinbc_a, const MF* robinbc_b, + const MF* robinbc_f) +{ + amrex::ignore_unused(amrlev, levelbcdata, robinbc_a, robinbc_b, robinbc_f); +} + +void MLCurlCurl::restriction (int amrlev, int cmglev, MF& crse, MF& fine) const +{ + IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[cmglev-1]; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + amrex::average_down_edges(fine[idim], crse[idim], ratio); + } +#if (AMREX_SPACEDIM == 2) + amrex::average_down_nodal(fine[2], crse[2], ratio); +#endif +} + +void MLCurlCurl::interpolation (int amrlev, int fmglev, MF& fine, + const MF& crse) const +{ + IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[fmglev]; + AMREX_ALWAYS_ASSERT(ratio == 2); + + for (int idim = 0; idim < 3; ++idim) { + auto const& finema = fine[idim].arrays(); + auto const& crsema = crse[idim].const_arrays(); + ParallelFor(fine[idim], [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k) + { + mlcurlcurl_interpadd(idim,i,j,k,finema[bno],crsema[bno]); + }); + } + Gpu::streamSynchronize(); +} + +void +MLCurlCurl::interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const +{ + amrex::ignore_unused(amrlev, fmglev, fine, crse); + amrex::Abort("MLCurlCurl::interpAssign: TODO"); +} + +void MLCurlCurl::interpolationAmr (int famrlev, MF& fine, const MF& crse, + IntVect const& nghost) const +{ + amrex::ignore_unused(famrlev, fine, crse, nghost); + amrex::Abort("MLCurlCurl::interpolationAmr: TODO"); +} + +void MLCurlCurl::averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs, + const MF& fine_sol, const MF& fine_rhs) +{ + amrex::ignore_unused(camrlev, crse_sol, crse_rhs, fine_sol, fine_rhs); + amrex::Abort("MLCurlCurl::averageDownSolutionRHS: TODO"); +} + +void +MLCurlCurl::apply (int amrlev, int mglev, MF& out, MF& in, BCMode /*bc_mode*/, + StateMode /*s_mode*/, const MLMGBndryT* /*bndry*/) const +{ + applyBC(amrlev, mglev, in); + + auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto const a = m_alpha; + auto const b = m_beta; + + int const dirichlet_xlo = getDirichlet(amrlev, mglev, 0, 0); + int const dirichlet_xhi = getDirichlet(amrlev, mglev, 0, 1); + int const dirichlet_ylo = getDirichlet(amrlev, mglev, 1, 0); + int const dirichlet_yhi = getDirichlet(amrlev, mglev, 1, 1); + int const dirichlet_zlo = getDirichlet(amrlev, mglev, 2, 0); + int const dirichlet_zhi = getDirichlet(amrlev, mglev, 2, 1); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(out[0],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& xbx = mfi.tilebox(out[0].ixType().toIntVect()); + Box const& ybx = mfi.tilebox(out[1].ixType().toIntVect()); + Box const& zbx = mfi.tilebox(out[2].ixType().toIntVect()); + auto const& xout = out[0].array(mfi); + auto const& yout = out[1].array(mfi); + auto const& zout = out[2].array(mfi); + auto const& xin = in[0].array(mfi); + auto const& yin = in[1].array(mfi); + auto const& zin = in[2].array(mfi); + amrex::ParallelFor(xbx, ybx, zbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (j == dirichlet_ylo || j == dirichlet_yhi || + k == dirichlet_zlo || k == dirichlet_zhi) { + xout(i,j,k) = Real(0.0); + } else { + mlcurlcurl_adotx_x(i,j,k,xout,xin,yin,zin,a,b,dxinv); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (i == dirichlet_xlo || i == dirichlet_xhi || + k == dirichlet_zlo || k == dirichlet_zhi) { + yout(i,j,k) = Real(0.0); + } else { + mlcurlcurl_adotx_y(i,j,k,yout,xin,yin,zin,a,b,dxinv); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (i == dirichlet_xlo || i == dirichlet_xhi || + j == dirichlet_ylo || j == dirichlet_yhi) { + zout(i,j,k) = Real(0.0); + } else { + mlcurlcurl_adotx_z(i,j,k,zout,xin,yin,zin,a,b,dxinv); + } + }); + } +} + +void MLCurlCurl::smooth (int amrlev, int mglev, MF& sol, const MF& rhs, + bool skip_fillboundary) const +{ + if (!skip_fillboundary) { + applyBC(amrlev, mglev, sol); + } + + smooth(amrlev, mglev, sol, rhs[0], 0); // Ex red + applyBC(amrlev, mglev, sol[0]); + + smooth(amrlev, mglev, sol, rhs[1], 0); // Ey red + applyBC(amrlev, mglev, sol[1]); + + smooth(amrlev, mglev, sol, rhs[2], 0); // Ez red + applyBC(amrlev, mglev, sol[2]); + + smooth(amrlev, mglev, sol, rhs[0], 1); // Ex black + applyBC(amrlev, mglev, sol[0]); + + smooth(amrlev, mglev, sol, rhs[1], 1); // Ey black +#if (AMREX_SPACEDIM == 3) + applyBC(amrlev, mglev, sol[1]); +#endif + + smooth(amrlev, mglev, sol, rhs[2], 1); // Ez black + + for (int idim = 0; idim < 3; ++idim) { + amrex::OverrideSync(sol[idim], getDotMask(amrlev,mglev,idim), + this->m_geom[amrlev][mglev].periodicity()); + } +} + +void MLCurlCurl::smooth (int amrlev, int mglev, MF& sol, MultiFab const& rhs, + int redblack) const +{ + auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto const a = m_alpha; + auto const b = m_beta; + + int const dirichlet_xlo = getDirichlet(amrlev, mglev, 0, 0); + int const dirichlet_xhi = getDirichlet(amrlev, mglev, 0, 1); + int const dirichlet_ylo = getDirichlet(amrlev, mglev, 1, 0); + int const dirichlet_yhi = getDirichlet(amrlev, mglev, 1, 1); + int const dirichlet_zlo = getDirichlet(amrlev, mglev, 2, 0); + int const dirichlet_zhi = getDirichlet(amrlev, mglev, 2, 1); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.tilebox(); + auto const& rh = rhs.const_array(mfi); + if (rhs.ixType() == sol[0].ixType()) { + auto const& ex = sol[0].array(mfi); + auto const& ey = sol[1].const_array(mfi); + auto const& ez = sol[2].const_array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (j != dirichlet_ylo && j != dirichlet_yhi && + k != dirichlet_zlo && k != dirichlet_zhi) { + mlcurlcurl_gsrb_x(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack); + } + }); + } else if (rhs.ixType() == sol[1].ixType()) { + auto const& ex = sol[0].const_array(mfi); + auto const& ey = sol[1].array(mfi); + auto const& ez = sol[2].const_array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (i != dirichlet_xlo && i != dirichlet_xhi && + k != dirichlet_zlo && k != dirichlet_zhi) { + mlcurlcurl_gsrb_y(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack); + } + }); + } else { + auto const& ex = sol[0].const_array(mfi); + auto const& ey = sol[1].const_array(mfi); + auto const& ez = sol[2].array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (i != dirichlet_xlo && i != dirichlet_xhi && + j != dirichlet_ylo && j != dirichlet_yhi) { + mlcurlcurl_gsrb_z(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack); + } + }); + } + } +} + +void MLCurlCurl::solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, + const MF* /*crse_bcdata*/) +{ + BL_PROFILE("MLCurlCurl::solutionResidual()"); + const int mglev = 0; + apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Solution); + compresid(amrlev, mglev, resid, b); +} + +void MLCurlCurl::correctionResidual (int amrlev, int mglev, MF& resid, MF& x, + const MF& b, BCMode bc_mode, + const MF* crse_bcdata) +{ + AMREX_ALWAYS_ASSERT(bc_mode != BCMode::Inhomogeneous && crse_bcdata == nullptr); + apply(amrlev, mglev, resid, x, BCMode::Homogeneous, StateMode::Correction); + compresid(amrlev, mglev, resid, b); +} + +void MLCurlCurl::compresid (int amrlev, int mglev, MF& resid, MF const& b) const +{ + int const dirichlet_xlo = getDirichlet(amrlev, mglev, 0, 0); + int const dirichlet_xhi = getDirichlet(amrlev, mglev, 0, 1); + int const dirichlet_ylo = getDirichlet(amrlev, mglev, 1, 0); + int const dirichlet_yhi = getDirichlet(amrlev, mglev, 1, 1); + int const dirichlet_zlo = getDirichlet(amrlev, mglev, 2, 0); + int const dirichlet_zhi = getDirichlet(amrlev, mglev, 2, 1); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(resid[0],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& xbx = mfi.tilebox(resid[0].ixType().toIntVect()); + Box const& ybx = mfi.tilebox(resid[1].ixType().toIntVect()); + Box const& zbx = mfi.tilebox(resid[2].ixType().toIntVect()); + auto const& resx = resid[0].array(mfi); + auto const& resy = resid[1].array(mfi); + auto const& resz = resid[2].array(mfi); + auto const& bx = b[0].array(mfi); + auto const& by = b[1].array(mfi); + auto const& bz = b[2].array(mfi); + amrex::ParallelFor(xbx, ybx, zbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (j == dirichlet_ylo || j == dirichlet_yhi || + k == dirichlet_zlo || k == dirichlet_zhi) { + resx(i,j,k) = Real(0.0); + } else { + resx(i,j,k) = bx(i,j,k) - resx(i,j,k); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (i == dirichlet_xlo || i == dirichlet_xhi || + k == dirichlet_zlo || k == dirichlet_zhi) { + resy(i,j,k) = Real(0.0); + } else { + resy(i,j,k) = by(i,j,k) - resy(i,j,k); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (i == dirichlet_xlo || i == dirichlet_xhi || + j == dirichlet_ylo || j == dirichlet_yhi) { + resz(i,j,k) = Real(0.0); + } else { + resz(i,j,k) = bz(i,j,k) - resz(i,j,k); + } + }); + } +} + +void MLCurlCurl::reflux (int crse_amrlev, + MF& res, const MF& crse_sol, const MF& crse_rhs, + MF& fine_res, MF& fine_sol, const MF& fine_rhs) const +{ + amrex::ignore_unused(crse_amrlev, res, crse_sol, crse_rhs, fine_res, + fine_sol, fine_rhs); + amrex::Abort("MLCurlCurl::reflux: TODO"); +} + +void MLCurlCurl::compFlux (int amrlev, const Array& fluxes, + MF& sol, Location loc) const +{ + amrex::ignore_unused(amrlev, fluxes, sol, loc); + amrex::Abort("MLCurlCurl::compFlux: TODO"); +} + +void MLCurlCurl::compGrad (int amrlev, const Array& grad, + MF& sol, Location loc) const +{ + amrex::ignore_unused(amrlev, grad, sol, loc); + amrex::Abort("MLCurlCurl::compGrad: TODO"); +} + +void MLCurlCurl::prepareForSolve () +{ + // xxxxx TODO MLCurlCurl::prepareForSolve +} + +Real MLCurlCurl::xdoty (int amrlev, int mglev, const MF& x, const MF& y, + bool local) const +{ + auto result = Real(0.0); + for (int idim = 0; idim < 3; ++idim) { + auto rtmp = MultiFab::Dot(getDotMask(amrlev,mglev,idim), + x[idim], 0, y[idim], 0, 1, 0, false); + result += rtmp; + } + if (!local) { + ParallelAllReduce::Sum(result, ParallelContext::CommunicatorSub()); + } + return result; +} + +Real MLCurlCurl::normInf (int /*amrlev*/, MF const& mf, bool local) const +{ + return amrex::norminf(mf, 0, m_ncomp, IntVect(0), local); +} + +void MLCurlCurl::averageDownAndSync (Vector& sol) const +{ + BL_PROFILE("MLCurlCurl::averageDownAndSync()"); + AMREX_ALWAYS_ASSERT(sol.size() == 1); + const int amrlev = 0; + const int mglev = 0; + for (int idim = 0; idim < 3; ++idim) { + amrex::OverrideSync(sol[amrlev][idim], getDotMask(amrlev,mglev,idim), + this->m_geom[amrlev][mglev].periodicity()); + } +} + +void MLCurlCurl::avgDownResAmr (int clev, MF& cres, MF const& fres) const +{ + amrex::ignore_unused(clev, cres, fres); + amrex::Abort("MLCurlCurl::avgDownResAmr: TODO"); +} + +void MLCurlCurl::avgDownResMG (int clev, MF& cres, MF const& fres) const +{ + amrex::ignore_unused(clev, cres, fres); + amrex::Abort("MLCurlCurl::avgDownResMG: TODO"); +} + +void MLCurlCurl::make (Vector >& mf, IntVect const& ng) const +{ + MLLinOpT::make(mf, ng); +} + +Array +MLCurlCurl::make (int amrlev, int mglev, IntVect const& ng) const +{ + MF r; + for (int idim = 0; idim < 3; ++idim) { + r[idim].define(amrex::convert(this->m_grids[amrlev][mglev], m_etype[idim]), + this->m_dmap[amrlev][mglev], m_ncomp, ng, MFInfo(), + *(this->m_factory)[amrlev][mglev]); + } + return r; +} + +Array +MLCurlCurl::makeAlias (MF const& mf) const +{ + MF r; + for (int idim = 0; idim < 3; ++idim) { + r[idim] = MultiFab(mf[idim], amrex::make_alias, 0, mf[idim].nComp()); + } + return r; +} + +Array +MLCurlCurl::makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const +{ + BoxArray cba = this->m_grids[amrlev][mglev]; + IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[mglev]; + cba.coarsen(ratio); + + MF r; + for (int idim = 0; idim < 3; ++idim) { + r[idim].define(amrex::convert(cba, m_etype[idim]), + this->m_dmap[amrlev][mglev], m_ncomp, ng); + } + return r; +} + +Array +MLCurlCurl::makeCoarseAmr (int famrlev, IntVect const& ng) const +{ + BoxArray cba = this->m_grids[famrlev][0]; + IntVect ratio(this->AMRRefRatio(famrlev-1)); + cba.coarsen(ratio); + + MF r; + for (int idim = 0; idim < 3; ++idim) { + r[idim].define(amrex::convert(cba, m_etype[idim]), + this->m_dmap[famrlev][0], m_ncomp, ng); + } + return r; +} + +void MLCurlCurl::applyBC (int amrlev, int mglev, MF& in) const +{ + Vector mfs{in.data(),&(in[1]),&(in[2])}; + FillBoundary(mfs, this->m_geom[amrlev][mglev].periodicity()); + for (auto& mf : in) { + applyPhysBC(amrlev, mglev, mf); + } +} + +void MLCurlCurl::applyBC (int amrlev, int mglev, MultiFab& mf) const +{ + mf.FillBoundary(this->m_geom[amrlev][mglev].periodicity()); + applyPhysBC(amrlev, mglev, mf); +} + +#ifdef AMREX_USE_GPU +struct MLCurlCurlBCTag { + Array4 fab; + Box bx; + Orientation face; + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + Box const& box() const noexcept { return bx; } +}; +#endif + +void MLCurlCurl::applyPhysBC (int amrlev, int mglev, MultiFab& mf) const +{ + auto const idxtype = mf.ixType(); + Box const domain = amrex::convert(this->m_geom[amrlev][mglev].Domain(), idxtype); + + MFItInfo mfi_info{}; + +#ifdef AMREX_USE_GPU + Vector tags; + mfi_info.DisableDeviceSync(); +#endif + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(mf,mfi_info); mfi.isValid(); ++mfi) { + auto const& vbx = mfi.validbox(); + auto const& a = mf.array(mfi); + for (OrientationIter oit; oit; ++oit) { + Orientation const face = oit(); + int const idim = face.coordDir(); + bool is_symmetric = face.isLow() + ? m_lobc[0][idim] == LinOpBCType::symmetry + : m_hibc[0][idim] == LinOpBCType::symmetry; + if (domain[face] == vbx[face] && is_symmetric) { + Box b = vbx; + int shift = face.isLow() ? -1 : 1; + b.setRange(idim, domain[face] + shift, 1); +#ifdef AMREX_USE_GPU + tags.emplace_back(MLCurlCurlBCTag{a,b,face}); +#else + amrex::LoopOnCpu(b, [&] (int i, int j, int k) + { + mlcurlcurl_bc_symmetry(i, j, k, face, idxtype, a); + }); +#endif + } + } + } + +#ifdef AMREX_USE_GPU + ParallelFor(tags, + [=] AMREX_GPU_DEVICE (int i, int j, int k, MLCurlCurlBCTag const& tag) noexcept + { + mlcurlcurl_bc_symmetry(i, j, k, tag.face, idxtype, tag.fab); + }); +#endif +} + +iMultiFab const& MLCurlCurl::getDotMask (int amrlev, int mglev, int idim) const +{ + if (m_dotmask[amrlev][mglev][idim] == nullptr) { + MultiFab tmp(amrex::convert(this->m_grids[amrlev][mglev], m_etype[idim]), + this->m_dmap[amrlev][mglev], 1, 0, MFInfo().SetAlloc(false)); + m_dotmask[amrlev][mglev][idim] = + tmp.OwnerMask(this->m_geom[amrlev][mglev].periodicity()); + } + return *m_dotmask[amrlev][mglev][idim]; +} + +int MLCurlCurl::getDirichlet (int amrlev, int mglev, int idim, int face) const +{ +#if (AMREX_SPACEDIM == 2) + if (idim == 2) { + return std::numeric_limits::lowest(); + } +#endif + + if (face == 0) { + if (m_lobc[0][idim] == LinOpBCType::Dirichlet) { + return m_geom[amrlev][mglev].Domain().smallEnd(idim); + } else { + return std::numeric_limits::lowest(); + } + } else { + if (m_hibc[0][idim] == LinOpBCType::Dirichlet) { + return m_geom[amrlev][mglev].Domain().bigEnd(idim) + 1; + } else { + return std::numeric_limits::max(); + } + } +} + +} diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H new file mode 100644 index 00000000000..6f7432aae4f --- /dev/null +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H @@ -0,0 +1,362 @@ +#ifndef AMREX_ML_CURL_CURL_K_H_ +#define AMREX_ML_CURL_CURL_K_H_ +#include + +#include + +namespace amrex { + +/* Index types + * E_x : (0,1,1) + * E_y : (1,0,1) + * E_z : (1,1,0) + * (curl E)_x : (1,0,0) + * (curl E)_y : (0,1,0) + * (curl E)_z : (0,0,1) + */ + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_adotx_x (int i, int j, int k, Array4 const& Ax, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Real alpha, Real beta, + GpuArray const& dxinv) +{ +#if (AMREX_SPACEDIM == 2) + amrex::ignore_unused(ez); + Real dyy = dxinv[1] * dxinv[1]; + Real dxy = dxinv[0] * dxinv[1]; + Real ccex = ex(i ,j ,k ) * dyy * Real(2.0) + - dyy * (ex(i ,j-1,k ) + + ex(i ,j+1,k )) + + dxy * (ey(i ,j-1,k ) + - ey(i ,j ,k ) + - ey(i+1,j-1,k ) + + ey(i+1,j ,k )); +#else + Real dyy = dxinv[1] * dxinv[1]; + Real dzz = dxinv[2] * dxinv[2]; + Real dxy = dxinv[0] * dxinv[1]; + Real dxz = dxinv[0] * dxinv[2]; + Real ccex = ex(i ,j ,k ) * (dyy+dzz)*Real(2.0) + - dyy * (ex(i ,j-1,k ) + + ex(i ,j+1,k )) + - dzz * (ex(i ,j ,k+1) + + ex(i ,j ,k-1)) + + dxy * (ey(i ,j-1,k ) + - ey(i ,j ,k ) + - ey(i+1,j-1,k ) + + ey(i+1,j ,k )) + + dxz * (ez(i ,j ,k-1) + - ez(i ,j ,k ) + - ez(i+1,j ,k-1) + + ez(i+1,j ,k )); +#endif + Ax(i,j,k) = alpha * ccex + beta * ex(i,j,k); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_adotx_y (int i, int j, int k, Array4 const& Ay, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Real alpha, Real beta, + GpuArray const& dxinv) +{ +#if (AMREX_SPACEDIM == 2) + amrex::ignore_unused(ez); + Real dxx = dxinv[0] * dxinv[0]; + Real dxy = dxinv[0] * dxinv[1]; + Real ccey = ey(i ,j ,k ) * dxx * Real(2.0) + - dxx * (ey(i-1,j ,k ) + + ey(i+1,j ,k )) + + dxy * (ex(i-1,j ,k ) + - ex(i ,j ,k ) + - ex(i-1,j+1,k ) + + ex(i ,j+1,k )); +#else + Real dxx = dxinv[0] * dxinv[0]; + Real dzz = dxinv[2] * dxinv[2]; + Real dxy = dxinv[0] * dxinv[1]; + Real dyz = dxinv[1] * dxinv[2]; + Real ccey = ey(i ,j ,k ) * (dxx+dzz)*Real(2.0) + - dxx * (ey(i-1,j ,k ) + + ey(i+1,j ,k )) + - dzz * (ey(i ,j ,k-1) + + ey(i ,j ,k+1)) + + dxy * (ex(i-1,j ,k ) + - ex(i ,j ,k ) + - ex(i-1,j+1,k ) + + ex(i ,j+1,k )) + + dyz * (ez(i ,j ,k-1) + - ez(i ,j ,k ) + - ez(i ,j+1,k-1) + + ez(i ,j+1,k )); +#endif + Ay(i,j,k) = alpha * ccey + beta * ey(i,j,k); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_adotx_z (int i, int j, int k, Array4 const& Az, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Real alpha, Real beta, + GpuArray const& dxinv) +{ +#if (AMREX_SPACEDIM == 2) + amrex::ignore_unused(ex,ey); + Real dxx = dxinv[0] * dxinv[0]; + Real dyy = dxinv[1] * dxinv[1]; + Real ccez = ez(i ,j ,k ) * (dxx+dyy)*Real(2.0) + - dxx * (ez(i-1,j ,k ) + + ez(i+1,j ,k )) + - dyy * (ez(i ,j-1,k ) + + ez(i ,j+1,k )); +#else + Real dxx = dxinv[0] * dxinv[0]; + Real dyy = dxinv[1] * dxinv[1]; + Real dxz = dxinv[0] * dxinv[2]; + Real dyz = dxinv[1] * dxinv[2]; + Real ccez = ez(i ,j ,k ) * (dxx+dyy)*Real(2.0) + - dxx * (ez(i-1,j ,k ) + + ez(i+1,j ,k )) + - dyy * (ez(i ,j-1,k ) + + ez(i ,j+1,k )) + + dxz * (ex(i-1,j ,k ) + - ex(i ,j ,k ) + - ex(i-1,j ,k+1) + + ex(i ,j ,k+1)) + + dyz * (ey(i ,j-1,k ) + - ey(i ,j ,k ) + - ey(i ,j-1,k+1) + + ey(i ,j ,k+1)); +#endif + Az(i,j,k) = alpha * ccez + beta * ez(i,j,k); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_gsrb_x (int i, int j, int k, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Array4 const& rhs, + Real alpha, Real beta, + GpuArray const& dxinv, int redblack) +{ + constexpr Real omega = Real(1.15); + + if ((i+j+k+redblack) % 2 != 0) { return; } + +#if (AMREX_SPACEDIM == 2) + amrex::ignore_unused(ez); + Real dyy = dxinv[1] * dxinv[1]; + Real dxy = dxinv[0] * dxinv[1]; + Real gamma = alpha * (dyy)*Real(2.0) + beta; + Real ccex = + - dyy * (ex(i ,j-1,k ) + + ex(i ,j+1,k )) + + dxy * (ey(i ,j-1,k ) + - ey(i ,j ,k ) + - ey(i+1,j-1,k ) + + ey(i+1,j ,k )); +#else + Real dyy = dxinv[1] * dxinv[1]; + Real dzz = dxinv[2] * dxinv[2]; + Real dxy = dxinv[0] * dxinv[1]; + Real dxz = dxinv[0] * dxinv[2]; + Real gamma = alpha * (dyy+dzz)*Real(2.0) + beta; + Real ccex = + - dyy * (ex(i ,j-1,k ) + + ex(i ,j+1,k )) + - dzz * (ex(i ,j ,k+1) + + ex(i ,j ,k-1)) + + dxy * (ey(i ,j-1,k ) + - ey(i ,j ,k ) + - ey(i+1,j-1,k ) + + ey(i+1,j ,k )) + + dxz * (ez(i ,j ,k-1) + - ez(i ,j ,k ) + - ez(i+1,j ,k-1) + + ez(i+1,j ,k )); +#endif + Real res = rhs(i,j,k) - (gamma*ex(i,j,k) + alpha*ccex); + ex(i,j,k) += omega/gamma * res; +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_gsrb_y (int i, int j, int k, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Array4 const& rhs, + Real alpha, Real beta, + GpuArray const& dxinv, int redblack) +{ + constexpr Real omega = Real(1.15); + + if ((i+j+k+redblack) % 2 != 0) { return; } + +#if (AMREX_SPACEDIM == 2) + amrex::ignore_unused(ez); + Real dxx = dxinv[0] * dxinv[0]; + Real dxy = dxinv[0] * dxinv[1]; + Real gamma = alpha * (dxx)*Real(2.0) + beta; + Real ccey = + - dxx * (ey(i-1,j ,k ) + + ey(i+1,j ,k )) + + dxy * (ex(i-1,j ,k ) + - ex(i ,j ,k ) + - ex(i-1,j+1,k ) + + ex(i ,j+1,k )); +#else + Real dxx = dxinv[0] * dxinv[0]; + Real dzz = dxinv[2] * dxinv[2]; + Real dxy = dxinv[0] * dxinv[1]; + Real dyz = dxinv[1] * dxinv[2]; + Real gamma = alpha * (dxx+dzz)*Real(2.0) + beta; + Real ccey = + - dxx * (ey(i-1,j ,k ) + + ey(i+1,j ,k )) + - dzz * (ey(i ,j ,k-1) + + ey(i ,j ,k+1)) + + dxy * (ex(i-1,j ,k ) + - ex(i ,j ,k ) + - ex(i-1,j+1,k ) + + ex(i ,j+1,k )) + + dyz * (ez(i ,j ,k-1) + - ez(i ,j ,k ) + - ez(i ,j+1,k-1) + + ez(i ,j+1,k )); +#endif + Real res = rhs(i,j,k) - (gamma*ey(i,j,k) + alpha*ccey); + ey(i,j,k) += omega/gamma * res; +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_gsrb_z (int i, int j, int k, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Array4 const& rhs, + Real alpha, Real beta, + GpuArray const& dxinv, int redblack) +{ + constexpr Real omega = Real(1.15); + + if ((i+j+k+redblack) % 2 != 0) { return; } + +#if (AMREX_SPACEDIM == 2) + amrex::ignore_unused(ex,ey); + Real dxx = dxinv[0] * dxinv[0]; + Real dyy = dxinv[1] * dxinv[1]; + Real gamma = alpha * (dxx+dyy)*Real(2.0) + beta; + Real ccez = + - dxx * (ez(i-1,j ,k ) + + ez(i+1,j ,k )) + - dyy * (ez(i ,j-1,k ) + + ez(i ,j+1,k )); +#else + Real dxx = dxinv[0] * dxinv[0]; + Real dyy = dxinv[1] * dxinv[1]; + Real dxz = dxinv[0] * dxinv[2]; + Real dyz = dxinv[1] * dxinv[2]; + Real gamma = alpha * (dxx+dyy)*Real(2.0) + beta; + Real ccez = + - dxx * (ez(i-1,j ,k ) + + ez(i+1,j ,k )) + - dyy * (ez(i ,j-1,k ) + + ez(i ,j+1,k )) + + dxz * (ex(i-1,j ,k ) + - ex(i ,j ,k ) + - ex(i-1,j ,k+1) + + ex(i ,j ,k+1)) + + dyz * (ey(i ,j-1,k ) + - ey(i ,j ,k ) + - ey(i ,j-1,k+1) + + ey(i ,j ,k+1)); +#endif + Real res = rhs(i,j,k) - (gamma*ez(i,j,k) + alpha*ccez); + ez(i,j,k) += omega/gamma * res; +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_interpadd (int dir, int i, int j, int k, + Array4 const& fine, + Array4 const& crse) +{ + int ic = amrex::coarsen(i,2); + int jc = amrex::coarsen(j,2); + int kc = amrex::coarsen(k,2); + if (dir == 0) { + bool j_is_odd = (jc*2 != j); + bool k_is_odd = (kc*2 != k); + if (j_is_odd && k_is_odd) { + fine(i,j,k) += Real(0.25) * + (crse(ic,jc,kc ) + crse(ic,jc+1,kc ) + + crse(ic,jc,kc+1) + crse(ic,jc+1,kc+1)); + } else if (j_is_odd) { + fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc+1,kc)); + } else if (k_is_odd) { + fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc,kc+1)); + } else { + fine(i,j,k) += crse(ic,jc,kc); + } + } else if (dir == 1) { + bool i_is_odd = (ic*2 != i); + bool k_is_odd = (kc*2 != k); + if (i_is_odd && k_is_odd) { + fine(i,j,k) += Real(0.25) * + (crse(ic ,jc,kc ) + crse(ic+1,jc,kc ) + + crse(ic ,jc,kc+1) + crse(ic+1,jc,kc+1)); + } else if (i_is_odd) { + fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic+1,jc,kc)); + } else if (k_is_odd) { + fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc,kc+1)); + } else { + fine(i,j,k) += crse(ic,jc,kc); + } + } else { + bool i_is_odd = (ic*2 != i); + bool j_is_odd = (jc*2 != j); + if (i_is_odd && j_is_odd) { + fine(i,j,k) += Real(0.25) * + (crse(ic ,jc ,kc) + crse(ic+1,jc ,kc) + + crse(ic ,jc+1,kc) + crse(ic+1,jc+1,kc)); + } else if (i_is_odd) { + fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic+1,jc,kc)); + } else if (j_is_odd) { + fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc+1,kc)); + } else { + fine(i,j,k) += crse(ic,jc,kc); + } + } +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_bc_symmetry (int i, int j, int k, Orientation face, IndexType it, + Array4 const& a) +{ + int const idir = face.coordDir(); + int offset = face.isLow() ? 1 : -1; + Real sign; + if (it.cellCentered(idir)) { + sign = Real(-1.0); + } else { + sign = Real(1.0); + offset *= 2; + } + + if (idir == 0) { + a(i,j,k) = sign * a(i+offset,j,k); + } else if (idir == 1) { + a(i,j,k) = sign * a(i,j+offset,k); + } else { + a(i,j,k) = sign * a(i,j,k+offset); + } +} + +} + +#endif diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H index 5fc8de10022..cf9eb34951a 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H @@ -17,6 +17,7 @@ #include #include +#include #include #include @@ -85,15 +86,6 @@ struct LinOpEnumType enum struct Location { FaceCenter, FaceCentroid, CellCenter, CellCentroid }; }; -template struct LinOpData {}; -// -template -struct LinOpData > > -{ - using fab_type = typename T::fab_type; - using value_type = typename T::value_type; -}; - template class MLMGT; template class MLCGSolverT; template class MLPoissonT; @@ -112,8 +104,8 @@ public: template friend class GMRESMLMGT; using MFType = MF; - using FAB = typename LinOpData::fab_type; - using RT = typename LinOpData::value_type; + using FAB = typename FabDataType::fab_type; + using RT = typename FabDataType::value_type; using BCType = LinOpBCType; using BCMode = LinOpEnumType::BCMode; @@ -625,6 +617,10 @@ protected: [[nodiscard]] bool isCellCentered () const noexcept { return m_ixtype == 0; } + [[nodiscard]] virtual IntVect getNGrowVectRestriction () const { + return isCellCentered() ? IntVect(0) : IntVect(1); + } + virtual void make (Vector >& mf, IntVect const& ng) const; [[nodiscard]] virtual MF make (int amrlev, int mglev, IntVect const& ng) const; @@ -1387,18 +1383,13 @@ template void MLLinOpT::make (Vector >& mf, IntVect const& ng) const { - if constexpr (IsMultiFabLike_v) { - mf.clear(); - mf.resize(m_num_amr_levels); - for (int alev = 0; alev < m_num_amr_levels; ++alev) { - mf[alev].resize(m_num_mg_levels[alev]); - for (int mlev = 0; mlev < m_num_mg_levels[alev]; ++mlev) { - mf[alev][mlev] = make(alev, mlev, ng); - } + mf.clear(); + mf.resize(m_num_amr_levels); + for (int alev = 0; alev < m_num_amr_levels; ++alev) { + mf[alev].resize(m_num_mg_levels[alev]); + for (int mlev = 0; mlev < m_num_mg_levels[alev]; ++mlev) { + mf[alev][mlev] = make(alev, mlev, ng); } - } else { - amrex::ignore_unused(mf, ng); - amrex::Abort("MLLinOpT::make: how did we get here?"); } } @@ -1486,8 +1477,8 @@ template void MLLinOpT::avgDownResMG (int clev, MF& cres, MF const& fres) const { - const int ncomp = this->getNComp(); if constexpr (amrex::IsFabArray::value) { + const int ncomp = this->getNComp(); #ifdef AMREX_USE_EB if (!fres.isAllRegular()) { if constexpr (std::is_same()) { diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMG.H b/Src/LinearSolvers/MLMG/AMReX_MLMG.H index 46069b7a26b..facc231ae86 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLMG.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLMG.H @@ -1020,7 +1020,7 @@ MLMGT::prepareForSolve (Vector const& a_sol, Vector const& makeSolvable(); } - IntVect ng = linop.isCellCentered() ? IntVect(0) : IntVect(1); + IntVect ng = linop.getNGrowVectRestriction(); if (cf_strategy == CFStrategy::ghostnodes) { ng = ng_rhs; } if (!solve_called) { linop.make(res, ng); diff --git a/Src/LinearSolvers/MLMG/Make.package b/Src/LinearSolvers/MLMG/Make.package index 769ccad4d26..a8f267d4c26 100644 --- a/Src/LinearSolvers/MLMG/Make.package +++ b/Src/LinearSolvers/MLMG/Make.package @@ -87,6 +87,12 @@ ifneq ($(BL_NO_FORT),TRUE) F90EXE_sources += AMReX_MLLinOp_nd.F90 endif +ifneq ($(DIM),1) + CEXE_headers += AMReX_MLCurlCurl.H + CEXE_sources += AMReX_MLCurlCurl.cpp + CEXE_headers += AMReX_MLCurlCurl_K.H +endif + VPATH_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/MLMG INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/LinearSolvers/MLMG diff --git a/Tests/LinearSolvers/ABecLaplacian_C/MyTestPlotfile.cpp b/Tests/LinearSolvers/ABecLaplacian_C/MyTestPlotfile.cpp index 6dbb2c55b59..6c48c90d247 100644 --- a/Tests/LinearSolvers/ABecLaplacian_C/MyTestPlotfile.cpp +++ b/Tests/LinearSolvers/ABecLaplacian_C/MyTestPlotfile.cpp @@ -78,7 +78,7 @@ MyTest::writePlotfile () const amrex::Print() << "Level " << ilev << " max-norm error: " << plotmf[ilev].norminf(3) << " 1-norm error: " << plotmf[ilev].norm1(3)*dvol - << " 2-norm error: " << plotmf[ilev].norm2(3)*dvol << std::endl; + << " 2-norm error: " << plotmf[ilev].norm2(3)*std::sqrt(dvol) << std::endl; } WriteMultiLevelPlotfile("plot", nlevels, amrex::GetVecOfConstPtrs(plotmf), diff --git a/Tests/LinearSolvers/CurlCurl/CMakeLists.txt b/Tests/LinearSolvers/CurlCurl/CMakeLists.txt new file mode 100644 index 00000000000..9dacdeb2fea --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/CMakeLists.txt @@ -0,0 +1,19 @@ +foreach(D IN LISTS AMReX_SPACEDIM) + if (D EQUAL 1) + return() + endif () + + set(_sources + main.cpp + MyTest.cpp + MyTest.H + initProb.cpp + initProb_K.H) + + set(_input_files inputs) + + setup_test(${D} _sources _input_files) + + unset(_sources) + unset(_input_files) +endforeach() diff --git a/Tests/LinearSolvers/CurlCurl/GNUmakefile b/Tests/LinearSolvers/CurlCurl/GNUmakefile new file mode 100644 index 00000000000..ab4cb8e762f --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/GNUmakefile @@ -0,0 +1,27 @@ +DEBUG = FALSE +USE_MPI = TRUE +USE_OMP = FALSE +COMP = gnu +DIM = 3 +BL_NO_FORT = TRUE + +USE_CUDA = FALSE +USE_SYCL = FALSE +USE_HIP = FALSE + +TINY_PROFILE = FALSE + +AMREX_HOME = ../../.. + +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + +include ./Make.package + +Pdirs := Base Boundary LinearSolvers/MLMG + +Ppack += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package) + +include $(Ppack) + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules + diff --git a/Tests/LinearSolvers/CurlCurl/Make.package b/Tests/LinearSolvers/CurlCurl/Make.package new file mode 100644 index 00000000000..8f306d8eb54 --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/Make.package @@ -0,0 +1,5 @@ + +CEXE_sources += main.cpp +CEXE_sources += MyTest.cpp initProb.cpp +CEXE_headers += MyTest.H +CEXE_headers += initProb_K.H diff --git a/Tests/LinearSolvers/CurlCurl/MyTest.H b/Tests/LinearSolvers/CurlCurl/MyTest.H new file mode 100644 index 00000000000..fc17a239ebb --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/MyTest.H @@ -0,0 +1,46 @@ +#ifndef MY_TEST_H_ +#define MY_TEST_H_ + +#include + +class MyTest +{ +public: + + MyTest (); + + void solve (); + +// public for cuda + void initProb (); + +private: + + void readParameters (); + void initData (); + + int n_cell = 128; + int max_grid_size = 64; + + // For MLMG solver + int verbose = 1; + int bottom_verbose = 0; + int max_iter = 100; + bool agglomeration = true; + bool consolidation = true; + int max_coarsening_level = 30; + + amrex::Geometry geom; + amrex::BoxArray grids; + amrex::DistributionMapping dmap; + + amrex::Array solution; + amrex::Array exact; + amrex::Array rhs; + + amrex::Real alpha_over_dx2 = 1.0; + amrex::Real alpha; + amrex::Real beta = 1.0; +}; + +#endif diff --git a/Tests/LinearSolvers/CurlCurl/MyTest.cpp b/Tests/LinearSolvers/CurlCurl/MyTest.cpp new file mode 100644 index 00000000000..fa72fcbca25 --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/MyTest.cpp @@ -0,0 +1,112 @@ +#include + +#include "MyTest.H" + +#include +#include + +using namespace amrex; + +MyTest::MyTest () +{ + readParameters(); + initData(); +} + +void +MyTest::solve () +{ + LPInfo info; + info.setAgglomeration(agglomeration); + info.setConsolidation(consolidation); + info.setMaxCoarseningLevel(max_coarsening_level); + + MLCurlCurl mlcc({geom}, {grids}, {dmap}, info); + + mlcc.setDomainBC({AMREX_D_DECL(LinOpBCType::symmetry, + LinOpBCType::Dirichlet, + LinOpBCType::Periodic)}, + {AMREX_D_DECL(LinOpBCType::Dirichlet, + LinOpBCType::symmetry, + LinOpBCType::Periodic)}); + + mlcc.setScalars(alpha, beta); + + MLMGT > mlmg(mlcc); + mlmg.setMaxIter(max_iter); + mlmg.setVerbose(verbose); + mlmg.setBottomVerbose(bottom_verbose); + for (auto& mf : solution) { + mf.setVal(Real(0)); + } + mlmg.solve({&solution}, {&rhs}, Real(1.0e-10), Real(0)); + + amrex::Print() << " Number of cells: " << n_cell << std::endl; + auto dvol = AMREX_D_TERM(geom.CellSize(0), *geom.CellSize(1), *geom.CellSize(2)); + Array names{"Ex", "Ey", "Ez"}; + for (int idim = 0; idim < 3; ++idim) { + MultiFab::Subtract(solution[idim], exact[idim], 0, 0, 1, 0); + auto e0 = solution[idim].norminf(); + auto e1 = solution[idim].norm1(0,geom.periodicity()); + e1 *= dvol; + auto e2 = solution[idim].norm2(0,geom.periodicity()); + e2 *= std::sqrt(dvol); + amrex::Print() << " " << names[idim] << " errors (max, L1, L2): " + << e0 << " " << e1 << " " << e2 << std::endl; + } +} + +void +MyTest::readParameters () +{ + ParmParse pp; + pp.query("n_cell", n_cell); + pp.query("max_grid_size", max_grid_size); + + pp.query("verbose", verbose); + pp.query("bottom_verbose", bottom_verbose); + pp.query("max_iter", max_iter); + pp.query("agglomeration", agglomeration); + pp.query("consolidation", consolidation); + pp.query("max_coarsening_level", max_coarsening_level); + + pp.query("alpha_over_dx2", alpha_over_dx2); + pp.query("beta", beta); +} + +void +MyTest::initData () +{ + RealBox rb({AMREX_D_DECL(0.,0.,0.)}, {AMREX_D_DECL(1.,1.,1.)}); + Array is_periodic{AMREX_D_DECL(0,0,1)}; + Geometry::Setup(&rb, 0, is_periodic.data()); + Box domain(IntVect(0), IntVect(n_cell-1)); + geom.define(domain); + + const Real dx = geom.CellSize(0); + alpha = alpha_over_dx2 * dx*dx; + + grids.define(domain); + grids.maxSize(max_grid_size); + dmap.define(grids); + + for (int idim = 0; idim < 3; ++idim) { + IntVect itype(1); +#if (AMREX_SPACEDIM == 2) + if (idim < AMREX_SPACEDIM) +#endif + { + itype[idim] = 0; + } + BoxArray const& ba = amrex::convert(grids, itype); + solution[idim].define(ba,dmap,1,1); + exact [idim].define(ba,dmap,1,1); + rhs [idim].define(ba,dmap,1,0); + } + + initProb(); + + for (int idim = 0; idim < 3; ++idim) { + exact[idim].LocalCopy(solution[idim], 0, 0, 1, IntVect(1)); + } +} diff --git a/Tests/LinearSolvers/CurlCurl/initProb.cpp b/Tests/LinearSolvers/CurlCurl/initProb.cpp new file mode 100644 index 00000000000..3acf779bd05 --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/initProb.cpp @@ -0,0 +1,32 @@ + +#include "MyTest.H" +#include "initProb_K.H" + +using namespace amrex; + +void +MyTest::initProb () +{ + const auto prob_lo = geom.ProbLoArray(); + const auto dx = geom.CellSizeArray(); + const auto a = alpha; + const auto b = beta; +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(rhs[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& gbx = mfi.tilebox(IntVect(1),IntVect(1)); + GpuArray,3> rhsfab{rhs[0].array(mfi), + rhs[1].array(mfi), + rhs[2].array(mfi)}; + GpuArray,3> solfab{solution[0].array(mfi), + solution[1].array(mfi), + solution[2].array(mfi)}; + amrex::ParallelFor(gbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + actual_init_prob(i,j,k,rhsfab,solfab,prob_lo,dx,a,b); + }); + } +} diff --git a/Tests/LinearSolvers/CurlCurl/initProb_K.H b/Tests/LinearSolvers/CurlCurl/initProb_K.H new file mode 100644 index 00000000000..66b7d2916cd --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/initProb_K.H @@ -0,0 +1,107 @@ +#ifndef INIT_PROB_K_H_ +#define INIT_PROB_K_H_ + +#include + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void actual_init_prob (int i, int j, int k, + amrex::GpuArray,3> const& rhs, + amrex::GpuArray,3> const& sol, + amrex::GpuArray const& problo, + amrex::GpuArray const& dx, + amrex::Real alpha, amrex::Real beta) +{ + using namespace amrex; + + constexpr Real pi = amrex::Math::pi(); + + Real xnd = problo[0] + Real(i)*dx[0]; + Real ynd = problo[1] + Real(j)*dx[1]; + Real xcc = xnd + Real(0.5)*dx[0]; + Real ycc = ynd + Real(0.5)*dx[1]; +#if (AMREX_SPACEDIM == 3) + Real znd = problo[2] + Real(k)*dx[2]; + Real zcc = znd + Real(0.5)*dx[2]; +#endif + + if (sol[0].contains(i,j,k)) { + Real x = xcc; + Real y = ynd; + Real Ex = std::sin(pi*x) * std::sin(Real(2.5)*pi*y); +#if (AMREX_SPACEDIM == 3) + Real z = znd; + Ex *= std::sin(Real(2.0)*pi*z + Real(1./3.)*pi); +#endif + sol[0](i,j,k) = Ex; + } + + if (sol[1].contains(i,j,k)) { + Real x = xnd; + Real y = ycc; + Real Ey = std::cos(Real(2.5)*pi*x) * std::sin(Real(3.)*pi*y); +#if (AMREX_SPACEDIM == 3) + Real z = znd; + Ey *= std::sin(Real(4.)*pi*z + Real(0.25)*pi); +#endif + sol[1](i,j,k) = Ey; + } + + if (sol[2].contains(i,j,k)) { + Real x = xnd; + Real y = ynd; + Real Ez = std::cos(Real(3.5)*pi*x) * std::sin(Real(3.5)*pi*y); +#if (AMREX_SPACEDIM == 3) + Real z = zcc; + Ez *= std::sin(Real(4.)*pi*z + Real(1./6.)*pi); +#endif + sol[2](i,j,k) = Ez; + } + + if (rhs[0].contains(i,j,k)) { + Real x = xcc; + Real y = ynd; +#if (AMREX_SPACEDIM == 2) + Real cce = Real(-7.5)*pi*pi*std::sin(Real(2.5)*pi*x)*std::cos(Real(3.)*pi*y) + + Real(6.25)*pi*pi*std::sin(pi*x)*std::sin(Real(2.5)*pi*y); +#else + Real z = znd; + Real cce = Real(-7.5)*pi*pi*std::sin(Real(2.5)*pi*x)*std::cos(Real(3.)*pi*y)*std::sin(Real(4.)*pi*z+Real(0.25)*pi) + + Real(6.25)*pi*pi*std::sin(pi*x)*std::sin(Real(2.5)*pi*y)*std::sin(Real(2.)*pi*z+Real(1./3.)*pi) + - Real(14.)*pi*pi*std::sin(Real(3.5)*pi*x)*std::sin(Real(3.5)*pi*y)*std::cos(Real(4.)*pi*z+Real(1./6.)*pi) + + Real(4.)*pi*pi*std::sin(pi*x)*std::sin(Real(2.5)*pi*y)*std::sin(Real(2.)*pi*z+Real(1./3.)*pi); +#endif + rhs[0](i,j,k) = alpha*cce + beta*sol[0](i,j,k); + } + + if (rhs[1].contains(i,j,k)) { + Real x = xnd; + Real y = ycc; +#if (AMREX_SPACEDIM == 2) + Real cce = Real(6.25)*pi*pi*std::cos(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y) + + Real(2.5)*pi*pi*std::cos(pi*x)*std::cos(Real(2.5)*pi*y); +#else + Real z = znd; + Real cce = Real(6.25)*pi*pi*std::cos(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y)*std::sin(Real(4.)*pi*z+Real(0.25)*pi) + + Real(2.5)*pi*pi*std::cos(pi*x)*std::cos(Real(2.5)*pi*y)*std::sin(Real(2.)*pi*z+Real(1./3.)*pi) + + Real(14.)*pi*pi*std::cos(Real(3.5)*pi*x)*std::cos(Real(3.5)*pi*y)*std::cos(Real(4.)*pi*z+Real(1./6.)*pi) + + Real(16.)*pi*pi*std::cos(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y)*std::sin(Real(4.)*pi*z+Real(0.25)*pi); +#endif + rhs[1](i,j,k) = alpha*cce + beta*sol[1](i,j,k); + } + + if (rhs[2].contains(i,j,k)) { + Real x = xnd; + Real y = ynd; +#if (AMREX_SPACEDIM == 2) + Real cce = Real(24.5)*pi*pi*std::cos(Real(3.5)*pi*x)*std::sin(Real(3.5)*pi*y); +#else + Real z = zcc; + Real cce = Real(24.5)*pi*pi*std::cos(Real(3.5)*pi*x)*std::sin(Real(3.5)*pi*y)*std::sin(Real(4.)*pi*z+Real(1./6.)*pi) + + Real(2.)*pi*pi*std::cos(pi*x)*std::sin(Real(2.5)*pi*y)*std::cos(Real(2.)*pi*z+Real(1./3.)*pi) + + Real(12.)*pi*pi*std::cos(Real(2.5)*pi*x)*std::cos(Real(3.)*pi*y)*std::cos(Real(4.)*pi*z+Real(0.25)*pi); +#endif + rhs[2](i,j,k) = alpha*cce + beta*sol[2](i,j,k); + } +} + +#endif diff --git a/Tests/LinearSolvers/CurlCurl/inputs b/Tests/LinearSolvers/CurlCurl/inputs new file mode 100644 index 00000000000..375562867ac --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/inputs @@ -0,0 +1,13 @@ + +n_cell = 128 +max_grid_size = 64 + +verbose = 2 +bottom_verbose = 2 + +alpha_over_dx2 = 1.0 +beta = 1.0 + +amrex.fpe_trap_invalid=1 +amrex.fpe_trap_zero=1 +amrex.fpe_trap_overflow=1 diff --git a/Tests/LinearSolvers/CurlCurl/main.cpp b/Tests/LinearSolvers/CurlCurl/main.cpp new file mode 100644 index 00000000000..769e28ce2b5 --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/main.cpp @@ -0,0 +1,17 @@ + +#include +#include +#include "MyTest.H" + +int main (int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + + { + BL_PROFILE("main"); + MyTest mytest; + mytest.solve(); + } + + amrex::Finalize(); +}