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/MLMG/AMReX_MLCurlCurl.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H index a74924625ec..3ec9b5d8e9d 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H @@ -116,25 +116,17 @@ public: // 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; + 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_x (int amrlev, int mglev, MultiFab& mf) const; - void applyPhysBC_y (int amrlev, int mglev, MultiFab& mf) const; - void applyPhysBC_z (int amrlev, int mglev, MultiFab& mf) const; + void applyPhysBC (int amrlev, int mglev, MultiFab& mf) const; private: void applyBC (int amrlev, int mglev, MF& in) const; - void applyBC_x (int amrlev, int mglev, MultiFab& mf) const; - void applyBC_y (int amrlev, int mglev, MultiFab& mf) const; - void applyBC_z (int amrlev, int mglev, MultiFab& mf) const; + void applyBC (int amrlev, int mglev, MultiFab& mf) const; [[nodiscard]] iMultiFab const& getDotMask (int amrlev, int mglev, int idim) const; diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp index 4240c7298db..2c86947af6b 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp @@ -153,22 +153,22 @@ void MLCurlCurl::smooth (int amrlev, int mglev, MF& sol, const MF& rhs, applyBC(amrlev, mglev, sol); } - smooth_x(amrlev, mglev, sol, rhs[0], 0); - applyBC_x(amrlev, mglev, sol[0]); + smooth(amrlev, mglev, sol, rhs[0], 0); // Ex red + applyBC(amrlev, mglev, sol[0]); - smooth_y(amrlev, mglev, sol, rhs[1], 0); - applyBC_y(amrlev, mglev, sol[1]); + smooth(amrlev, mglev, sol, rhs[1], 0); // Ey red + applyBC(amrlev, mglev, sol[1]); - smooth_z(amrlev, mglev, sol, rhs[2], 0); - applyBC_z(amrlev, mglev, sol[2]); + smooth(amrlev, mglev, sol, rhs[2], 0); // Ez red + applyBC(amrlev, mglev, sol[2]); - smooth_x(amrlev, mglev, sol, rhs[0], 1); - applyBC_x(amrlev, mglev, sol[0]); + smooth(amrlev, mglev, sol, rhs[0], 1); // Ex black + applyBC(amrlev, mglev, sol[0]); - smooth_y(amrlev, mglev, sol, rhs[1], 1); - applyBC_y(amrlev, mglev, sol[1]); + smooth(amrlev, mglev, sol, rhs[1], 1); // Ey black + applyBC(amrlev, mglev, sol[1]); - smooth_z(amrlev, mglev, sol, rhs[2], 1); + smooth(amrlev, mglev, sol, rhs[2], 1); // Ez black for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { amrex::OverrideSync(sol[idim], getDotMask(amrlev,mglev,idim), @@ -176,79 +176,9 @@ void MLCurlCurl::smooth (int amrlev, int mglev, MF& sol, const MF& rhs, } } -void MLCurlCurl::smooth_x (int amrlev, int mglev, MF& sol, MultiFab const& rhs, - int redblack) const +void MLCurlCurl::smooth (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; - - 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(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) - { - 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); - } - }); - } -} - -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; - - int const dirichlet_xlo = getDirichlet(amrlev, mglev, 0, 0); - int const dirichlet_xhi = getDirichlet(amrlev, mglev, 0, 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(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) - { - 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); - } - }); - } -} - -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; @@ -257,24 +187,50 @@ void MLCurlCurl::smooth_z (int amrlev, int mglev, MF& sol, MultiFab const& rhs, 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(sol[2],TilingIfNotGPU()); mfi.isValid(); ++mfi) + for (MFIter mfi(rhs,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) - { - 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); - } - }); + 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); + } + }); + } } } @@ -483,42 +439,45 @@ void MLCurlCurl::applyBC (int amrlev, int mglev, MF& in) const { Vector mfs{AMREX_D_DECL(in.data(),&(in[1]),&(in[2]))}; FillBoundary(mfs, this->m_geom[amrlev][mglev].periodicity()); - applyPhysBC_x(amrlev, mglev, in[0]); - applyPhysBC_y(amrlev, mglev, in[1]); - applyPhysBC_z(amrlev, mglev, in[2]); + applyPhysBC(amrlev, mglev, in[0]); + applyPhysBC(amrlev, mglev, in[1]); + applyPhysBC(amrlev, mglev, in[2]); } -void MLCurlCurl::applyBC_x (int amrlev, int mglev, MultiFab& mf) const +void MLCurlCurl::applyBC (int amrlev, int mglev, MultiFab& mf) const { mf.FillBoundary(this->m_geom[amrlev][mglev].periodicity()); - applyPhysBC_x(amrlev, mglev, mf); -} - -void MLCurlCurl::applyBC_y (int amrlev, int mglev, MultiFab& mf) const -{ - mf.FillBoundary(this->m_geom[amrlev][mglev].periodicity()); - applyPhysBC_y(amrlev, mglev, mf); -} - -void MLCurlCurl::applyBC_z (int amrlev, int mglev, MultiFab& mf) const -{ - mf.FillBoundary(this->m_geom[amrlev][mglev].periodicity()); - applyPhysBC_z(amrlev, mglev, mf); -} - -void MLCurlCurl::applyPhysBC_x (int amrlev, int mglev, MultiFab& mf) const -{ - -} - -void MLCurlCurl::applyPhysBC_y (int amrlev, int mglev, MultiFab& mf) const -{ - -} - -void MLCurlCurl::applyPhysBC_z (int amrlev, int mglev, MultiFab& mf) const -{ - + applyPhysBC(amrlev, mglev, mf); +} + +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); + for (MFIter mfi(mf); 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 + static_assert(false, "MLCurlCurl: todo"); +#else + amrex::LoopOnCpu(b, [&] (int i, int j, int k) + { + mlcurlcurl_bc_symmetry(i, j, k, face, idxtype, a); + }); +#endif + } + } + } } iMultiFab const& MLCurlCurl::getDotMask (int amrlev, int mglev, int idim) const diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H index 72965e0f3fa..4ff3e561a5d 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H @@ -257,6 +257,29 @@ void mlcurlcurl_interpadd (int dir, int i, int j, int k, } } +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()) { + 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/Tests/LinearSolvers/CurlCurl/MyTest.H b/Tests/LinearSolvers/CurlCurl/MyTest.H index 5d81f5e4026..68c8bca0058 100644 --- a/Tests/LinearSolvers/CurlCurl/MyTest.H +++ b/Tests/LinearSolvers/CurlCurl/MyTest.H @@ -38,7 +38,7 @@ private: amrex::Array exact; amrex::Array rhs; - amrex::Real alpha_over_dx2 = 1.e-1; + amrex::Real alpha_over_dx2 = 1.0; amrex::Real alpha; amrex::Real beta = 1.0; }; diff --git a/Tests/LinearSolvers/CurlCurl/MyTest.cpp b/Tests/LinearSolvers/CurlCurl/MyTest.cpp index 996faab5ee8..5ce1a483070 100644 --- a/Tests/LinearSolvers/CurlCurl/MyTest.cpp +++ b/Tests/LinearSolvers/CurlCurl/MyTest.cpp @@ -23,12 +23,12 @@ MyTest::solve () MLCurlCurl mlcc({geom}, {grids}, {dmap}, info); - mlcc.setDomainBC({AMREX_D_DECL(LinOpBCType::Dirichlet, + mlcc.setDomainBC({AMREX_D_DECL(LinOpBCType::symmetry, LinOpBCType::Periodic, LinOpBCType::Dirichlet)}, {AMREX_D_DECL(LinOpBCType::Dirichlet, LinOpBCType::Periodic, - LinOpBCType::Dirichlet)}); + LinOpBCType::symmetry)}); mlcc.setScalars(alpha, beta);