From 77ec6dbb074cd64c87000c24bfcdfe105e682cde 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 3D, where E is an array of MultiFabs on edges, and alpha and beta are scalars. Currently this is just a prototype. It only supports periodic boundaries. --- Src/Base/AMReX_FabDataType.H | 27 ++ Src/Base/CMakeLists.txt | 1 + Src/Base/Make.package | 1 + Src/Boundary/AMReX_FabSet.H | 5 +- Src/LinearSolvers/CMakeLists.txt | 9 + Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H | 142 +++++++ Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp | 401 ++++++++++++++++++++ Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H | 262 +++++++++++++ Src/LinearSolvers/MLMG/AMReX_MLLinOp.H | 33 +- Src/LinearSolvers/MLMG/Make.package | 6 + 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 | 96 +++++ Tests/LinearSolvers/CurlCurl/initProb.cpp | 34 ++ Tests/LinearSolvers/CurlCurl/initProb_K.H | 75 ++++ Tests/LinearSolvers/CurlCurl/inputs | 13 + Tests/LinearSolvers/CurlCurl/main.cpp | 17 + 19 files changed, 1194 insertions(+), 25 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/LinearSolvers/CMakeLists.txt b/Src/LinearSolvers/CMakeLists.txt index c2851d49959..94419bb201a 100644 --- a/Src/LinearSolvers/CMakeLists.txt +++ b/Src/LinearSolvers/CMakeLists.txt @@ -64,6 +64,15 @@ foreach(D IN LISTS AMReX_SPACEDIM) ) endif () + if (D EQUAL 3) + 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..ab70b0f1049 --- /dev/null +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H @@ -0,0 +1,142 @@ +#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 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: Make sure isCellCentered in MLMG makes sense. + */ +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; + + 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_x (int amrlev, int mglev, MF& sol, MultiFab const& rhs, + int redblack) const; + void smooth_y (int amrlev, int mglev, MF& sol, MultiFab const& rhs, + int redblack) const; + void smooth_z (int amrlev, int mglev, MF& sol, MultiFab const& rhs, + int redblack) const; + +private: + + void applyBC (int amrlev, int mglev, MF& in, BCMode bc_mode, StateMode s_mode) const; + iMultiFab const& getDotMask (int amrlev, int mglev, int idim) const; + + RT m_alpha = std::numeric_limits::lowest(); + RT m_beta = std::numeric_limits::lowest(); + Array m_etype{IntVect(0,1,1), + IntVect(1,0,1), + IntVect(1,1,0)}; + mutable Vector,AMREX_SPACEDIM>>> 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..63eb09d54d2 --- /dev/null +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp @@ -0,0 +1,401 @@ + +#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); + } +} + +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 < AMREX_SPACEDIM; ++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, bc_mode, s_mode); + + auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto const a = m_alpha; + auto const b = m_beta; + +#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) + { + mlcurlcurl_adotx_x(i,j,k,xout,xin,yin,zin,a,b,dxinv); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + mlcurlcurl_adotx_y(i,j,k,yout,xin,yin,zin,a,b,dxinv); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + 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 +{ + // xxxxx TODO: note that applyBC unnecessarily fill all three multifabs + + if (!skip_fillboundary) { + applyBC(amrlev, mglev, sol, BCMode::Homogeneous, StateMode::Correction); + } + smooth_x(amrlev, mglev, sol, rhs[0], 0); + + applyBC(amrlev, mglev, sol, BCMode::Homogeneous, StateMode::Correction); + smooth_y(amrlev, mglev, sol, rhs[1], 0); + + applyBC(amrlev, mglev, sol, BCMode::Homogeneous, StateMode::Correction); + smooth_z(amrlev, mglev, sol, rhs[2], 0); + + applyBC(amrlev, mglev, sol, BCMode::Homogeneous, StateMode::Correction); + smooth_x(amrlev, mglev, sol, rhs[0], 1); + + applyBC(amrlev, mglev, sol, BCMode::Homogeneous, StateMode::Correction); + smooth_y(amrlev, mglev, sol, rhs[1], 1); + + applyBC(amrlev, mglev, sol, BCMode::Homogeneous, StateMode::Correction); + smooth_z(amrlev, mglev, sol, rhs[2], 1); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + amrex::OverrideSync(sol[idim], getDotMask(amrlev,mglev,idim), + this->m_geom[amrlev][mglev].periodicity()); + } +} + +void MLCurlCurl::smooth_x (int amrlev, int mglev, MF& sol, MultiFab const& rhs, + int redblack) const +{ + amrex::ignore_unused(amrlev,mglev); + + auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto const a = m_alpha; + auto const b = m_beta; + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(sol[0],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.tilebox(); + auto const& ex = sol[0].array(mfi); + auto const& ey = sol[1].const_array(mfi); + auto const& ez = sol[2].const_array(mfi); + auto const& rh = rhs.const_array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + mlcurlcurl_gsrb_x(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack); + }); + } +} + +void MLCurlCurl::smooth_y (int amrlev, int mglev, MF& sol, MultiFab const& rhs, + int redblack) const +{ + amrex::ignore_unused(amrlev,mglev); + + auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto const a = m_alpha; + auto const b = m_beta; + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(sol[1],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.tilebox(); + auto const& ex = sol[0].const_array(mfi); + auto const& ey = sol[1]. array(mfi); + auto const& ez = sol[2].const_array(mfi); + auto const& rh = rhs.const_array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + mlcurlcurl_gsrb_y(i,j,k,ex,ey,ez,rh,a,b,dxinv,redblack); + }); + } +} + +void MLCurlCurl::smooth_z (int amrlev, int mglev, MF& sol, MultiFab const& rhs, + int redblack) const +{ + amrex::ignore_unused(amrlev,mglev); + + auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto const a = m_alpha; + auto const b = m_beta; + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(sol[2],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box const& bx = mfi.tilebox(); + auto const& ex = sol[0].const_array(mfi); + auto const& ey = sol[1].const_array(mfi); + auto const& ez = sol[2]. array(mfi); + auto const& rh = rhs.const_array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + 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); + amrex::Xpay(resid, RT(-1.0), b, 0, 0, m_ncomp, IntVect(0)); +} + +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); + amrex::Xpay(resid, RT(-1.0), b, 0, 0, m_ncomp, IntVect(0)); +} + +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 < AMREX_SPACEDIM; ++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 < AMREX_SPACEDIM; ++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 < AMREX_SPACEDIM; ++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 < AMREX_SPACEDIM; ++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 < AMREX_SPACEDIM; ++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 < AMREX_SPACEDIM; ++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, BCMode /*bc_mode*/, + StateMode /*s_mode*/) const +{ + Vector mfs{AMREX_D_DECL(in.data(),&(in[1]),&(in[2]))}; + FillBoundary(mfs, this->m_geom[amrlev][mglev].periodicity()); +} + +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]; +} + +} diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H new file mode 100644 index 00000000000..72965e0f3fa --- /dev/null +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H @@ -0,0 +1,262 @@ +#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) +{ + 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 )); + 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) +{ + 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 )); + 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) +{ + 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)); + 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; } + + 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 )); + 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; } + + 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 )); + 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; } + + 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)); + 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); + } + } +} + +} + +#endif diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H index f0dca07f3ab..52b0be4937a 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; @@ -109,8 +101,8 @@ public: template friend class MLPoissonT; template friend class MLABecLaplacianT; - 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; @@ -1384,18 +1376,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?"); } } @@ -1483,8 +1470,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/Make.package b/Src/LinearSolvers/MLMG/Make.package index d66d64ec0eb..03146d64808 100644 --- a/Src/LinearSolvers/MLMG/Make.package +++ b/Src/LinearSolvers/MLMG/Make.package @@ -85,5 +85,11 @@ ifneq ($(BL_NO_FORT),TRUE) F90EXE_sources += AMReX_MLLinOp_nd.F90 endif +ifeq ($(DIM),3) + 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/CurlCurl/CMakeLists.txt b/Tests/LinearSolvers/CurlCurl/CMakeLists.txt new file mode 100644 index 00000000000..afc06542acd --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/CMakeLists.txt @@ -0,0 +1,19 @@ +foreach(D IN LISTS AMReX_SPACEDIM) + if (NOT (D EQUAL 3)) + 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..3e52a8d11f0 --- /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: // make it 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.e-1; + 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..3bba1ccb058 --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/MyTest.cpp @@ -0,0 +1,96 @@ +#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::Periodic, + LinOpBCType::Periodic, + LinOpBCType::Periodic)}, + {AMREX_D_DECL(LinOpBCType::Periodic, + LinOpBCType::Periodic, + LinOpBCType::Periodic)}); + + + mlcc.setLevelBC(0, &exact); + + 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)); +} + +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(1,1,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 < AMREX_SPACEDIM; ++idim) { + IntVect itype(1); + 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 < AMREX_SPACEDIM; ++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..0d8d98d36cd --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/initProb.cpp @@ -0,0 +1,34 @@ + +#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,AMREX_SPACEDIM> rhsfab + {AMREX_D_DECL(rhs[0].array(mfi), + rhs[1].array(mfi), + rhs[2].array(mfi))}; + GpuArray,AMREX_SPACEDIM> solfab + {AMREX_D_DECL(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..0401b5c9b7a --- /dev/null +++ b/Tests/LinearSolvers/CurlCurl/initProb_K.H @@ -0,0 +1,75 @@ +#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,AMREX_SPACEDIM> const& rhs, + amrex::GpuArray,AMREX_SPACEDIM> const& sol, + amrex::GpuArray const& problo, + amrex::GpuArray const& dx, + amrex::Real alpha, amrex::Real beta) +{ + using namespace amrex; + + constexpr amrex::Real tpi = Real(2.)*amrex::Math::pi(); + constexpr amrex::Real fpi = Real(4.)*amrex::Math::pi(); + + Real xnd = problo[0] + Real(i)*dx[0]; + Real ynd = problo[1] + Real(j)*dx[1]; + Real znd = problo[2] + Real(k)*dx[2]; + Real xcc = xnd + Real(0.5)*dx[0]; + Real ycc = ynd + Real(0.5)*dx[1]; + Real zcc = znd + Real(0.5)*dx[2]; + + if (sol[0].contains(i,j,k)) { + Real x = xcc; + Real y = ynd; + Real z = znd; + sol[0](i,j,k) = std::cos(tpi*x) * std::sin(fpi*y) * std::sin(fpi*z); + } + + if (sol[1].contains(i,j,k)) { + Real x = xnd; + Real y = ycc; + Real z = znd; + sol[1](i,j,k) = std::sin(fpi*x) * std::cos(tpi*y) * std::sin(fpi*z); + } + + if (sol[2].contains(i,j,k)) { + Real x = xnd; + Real y = ynd; + Real z = zcc; + sol[2](i,j,k) = std::sin(fpi*x) * std::sin(fpi*y) * std::cos(tpi*z); + } + + if (rhs[0].contains(i,j,k)) { + Real x = xcc; + Real y = ynd; + Real z = znd; + rhs[0](i,j,k) = (beta + alpha*Real(4.0)*tpi*fpi) * sol[0](i,j,k) + - alpha*tpi*fpi * std::cos(fpi*x) * (std::sin(tpi*y) * std::sin(fpi*z) + + std::sin(fpi*y) * std::sin(tpi*z)); + } + + if (rhs[1].contains(i,j,k)) { + Real x = xnd; + Real y = ycc; + Real z = znd; + rhs[1](i,j,k) = (beta + alpha*Real(4.0)*tpi*fpi) * sol[1](i,j,k) + - alpha*tpi*fpi * std::cos(fpi*y) * (std::sin(tpi*x) * std::sin(fpi*z) + + std::sin(fpi*x) * std::sin(tpi*z)); + } + + if (rhs[2].contains(i,j,k)) { + Real x = xnd; + Real y = ynd; + Real z = zcc; + rhs[2](i,j,k) = (beta + alpha*Real(4.0)*tpi*fpi) * sol[2](i,j,k) + - alpha*tpi*fpi * std::cos(fpi*z) * (std::sin(tpi*x) * std::sin(fpi*y) + + std::sin(fpi*x) * std::sin(tpi*y)); + } +} + +#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(); +}