From 3a8ce866a5e6c6cc9ba864751322e3857966bb0c Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 19 Dec 2023 12:46:12 -0800 Subject: [PATCH] New Linear Solver: Curl of Curl Add a new linear solver for curl (alpha curl E) + beta E = rhs in 2D & 3D, where E is an array of 3 MultiFabs on edges, and alpha and beta are scalars. It supports periodic, homogeneous Dirichlet, and symmetry boundaries. At the symmetry boundary, the normal component changes the sign, whereas the transverse components do not. --- Src/Base/AMReX_FabDataType.H | 27 + Src/Base/CMakeLists.txt | 1 + Src/Base/Make.package | 1 + Src/Boundary/AMReX_FabSet.H | 5 +- Src/Boundary/AMReX_LO_BCTYPES.H | 3 +- Src/Boundary/AMReX_LO_BCTYPES.cpp | 5 + Src/LinearSolvers/CMakeLists.txt | 9 + Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H | 155 +++++ Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp | 552 ++++++++++++++++++ Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H | 362 ++++++++++++ Src/LinearSolvers/MLMG/AMReX_MLLinOp.H | 37 +- Src/LinearSolvers/MLMG/AMReX_MLMG.H | 2 +- Src/LinearSolvers/MLMG/Make.package | 6 + .../ABecLaplacian_C/MyTestPlotfile.cpp | 2 +- Tests/LinearSolvers/CurlCurl/CMakeLists.txt | 19 + Tests/LinearSolvers/CurlCurl/GNUmakefile | 27 + Tests/LinearSolvers/CurlCurl/Make.package | 5 + Tests/LinearSolvers/CurlCurl/MyTest.H | 46 ++ Tests/LinearSolvers/CurlCurl/MyTest.cpp | 112 ++++ Tests/LinearSolvers/CurlCurl/initProb.cpp | 32 + Tests/LinearSolvers/CurlCurl/initProb_K.H | 107 ++++ Tests/LinearSolvers/CurlCurl/inputs | 13 + Tests/LinearSolvers/CurlCurl/main.cpp | 17 + 23 files changed, 1517 insertions(+), 28 deletions(-) create mode 100644 Src/Base/AMReX_FabDataType.H create mode 100644 Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H create mode 100644 Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp create mode 100644 Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H create mode 100644 Tests/LinearSolvers/CurlCurl/CMakeLists.txt create mode 100644 Tests/LinearSolvers/CurlCurl/GNUmakefile create mode 100644 Tests/LinearSolvers/CurlCurl/Make.package create mode 100644 Tests/LinearSolvers/CurlCurl/MyTest.H create mode 100644 Tests/LinearSolvers/CurlCurl/MyTest.cpp create mode 100644 Tests/LinearSolvers/CurlCurl/initProb.cpp create mode 100644 Tests/LinearSolvers/CurlCurl/initProb_K.H create mode 100644 Tests/LinearSolvers/CurlCurl/inputs create mode 100644 Tests/LinearSolvers/CurlCurl/main.cpp 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(); +}