diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_1D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_1D_K.H index e74041edbc..4f982e0775 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_1D_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLap_1D_K.H @@ -5,164 +5,425 @@ namespace amrex { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_set_nodal_mask (int /*i*/, int /*j*/, int /*k*/, Array4 const& /*nmsk*/, - Array4 const& /*cmsk*/) noexcept -{} +void mlndlap_set_nodal_mask (int i, int, int, Array4 const& nmsk, + Array4 const& cmsk) noexcept +{ + using namespace nodelap_detail; + + int s = cmsk(i-1,0,0) + cmsk(i,0,0); + if (s == 2*crse_cell) { + nmsk(i,0,0) = crse_node; + } else if (s == 2*fine_cell) { + nmsk(i,0,0) = fine_node; + } else { + nmsk(i,0,0) = crse_fine_node; + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_set_dirichlet_mask (Box const&, Array4 const&, - Array4 const&, Box const&, - GpuArray const&, - GpuArray const&) noexcept -{} +void mlndlap_set_dirichlet_mask (Box const& bx, Array4 const& dmsk, + Array4 const& omsk, Box const& dom, + GpuArray const& bclo, + GpuArray const& bchi) noexcept +{ + const auto lo = bx.smallEnd(0); + const auto hi = bx.bigEnd(0); + AMREX_PRAGMA_SIMD + for (int i = lo; i <= hi; ++i) { + if (!dmsk(i,0,0)) { + dmsk(i,0,0) = (omsk(i-1,0,0) == 1 || omsk(i,0,0) == 1); + } + } + + const auto domlo = dom.smallEnd(0); + const auto domhi = dom.bigEnd(0); + + if (bclo[0] == LinOpBCType::Dirichlet && lo == domlo) { + dmsk(lo,0,0) = 1; + } + + if (bchi[0] == LinOpBCType::Dirichlet && hi == domhi) { + dmsk(hi,0,0) = 1; + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_set_dot_mask (Box const&, Array4 const&, - Array4 const&, Box const&, - GpuArray const&, - GpuArray const&) noexcept -{} +void mlndlap_set_dot_mask (Box const& bx, Array4 const& dmsk, + Array4 const& omsk, Box const& dom, + GpuArray const& bclo, + GpuArray const& bchi) noexcept +{ + const auto lo = bx.smallEnd(0); + const auto hi = bx.bigEnd(0); + + AMREX_PRAGMA_SIMD + for (int i = lo; i <= hi; ++i) { + dmsk(i,0,0) = static_cast(omsk(i,0,0)); + } + + const auto domlo = dom.smallEnd(0); + const auto domhi = dom.bigEnd(0); + + if ((bclo[0] == LinOpBCType::Neumann || bclo[0] == LinOpBCType::inflow) + && lo == domlo) + { + dmsk(lo,0,0) *= Real(0.5); + } + + if ((bchi[0] == LinOpBCType::Neumann || bchi[0] == LinOpBCType::inflow) + && hi == domhi) + { + dmsk(hi,0,0) *= Real(0.5); + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_zero_fine (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&, int) noexcept -{} +void mlndlap_zero_fine (int i, int, int, Array4 const& phi, + Array4 const& msk, int fine_flag) noexcept +{ + // Testing if the node is covered by a fine level in computing + // coarse sync residual + if (msk(i-1,0,0) == fine_flag && + msk(i ,0,0) == fine_flag) + { + phi(i,0,0) = Real(0.0); + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_avgdown_coeff_x (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&) noexcept -{} +void mlndlap_avgdown_coeff_x (int i, int, int, Array4 const& crse, + Array4 const& fine) noexcept +{ + Real a = fine(2*i ,0,0); + Real b = fine(2*i+1,0,0); + crse(i,0,0) = Real(2.0)*a*b/(a+b); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_semi_avgdown_coeff (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&, int) noexcept -{} +void mlndlap_semi_avgdown_coeff (int i, int j, int k, Array4 const& crse, + Array4 const& fine, int) noexcept +{ + mlndlap_avgdown_coeff_x(i,j,k,crse,fine); +} template -inline void mlndlap_bc_doit (Box const&, Array4 const&, Box const&, - GpuArray const&, - GpuArray const&) noexcept -{} +void mlndlap_bc_doit (Box const& vbx, Array4 const& a, Box const& domain, + GpuArray const& bflo, + GpuArray const& bfhi) noexcept +{ + Box gdomain = domain; + int const idim = 0; + if (! bflo[idim]) { gdomain.growLo(idim,1); } + if (! bfhi[idim]) { gdomain.growHi(idim,1); } + + if (gdomain.strictly_contains(vbx)) { return; } + + const int offset = domain.cellCentered() ? 0 : 1; + + const auto dlo = domain.smallEnd(0); + const auto dhi = domain.bigEnd(0); + + Box const& sbox = amrex::grow(vbx,1); + AMREX_HOST_DEVICE_FOR_3D(sbox, i, j, k, + { + if (! gdomain.contains(IntVect(i))) { + if (i == dlo-1 && bflo[0]) + { + a(i,0,0) = a(i+1+offset, j, k); + } + else if (i == dhi+1 && bfhi[0]) + { + a(i,0,0) = a(i-1-offset, j, k); + } + } + }); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -Real mlndlap_adotx_c (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Real, Array4 const&, - GpuArray const&) noexcept -{ return Real(0.0); } +Real mlndlap_adotx_c (int i, int, int, Array4 const& x, + Real sigma, Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + if (msk(i,0,0)) { + return Real(0.0); + } else { + Real y = x(i-1,0,0) + + x(i+1,0,0) + + x(i ,0,0)*Real(-2.0); + return dxinv[0]*dxinv[0] * y * sigma; + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -Real mlndlap_adotx_ha (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&, Array4 const&, - GpuArray const&) noexcept -{ return Real(0.0); } +Real mlndlap_adotx_ha (int i, int, int, Array4 const& x, + Array4 const& sx, Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + if (msk(i,0,0)) { + return Real(0.0); + } else { + Real y = (x(i+1,0,0) - x(i ,0,0)) * sx(i ,0,0) + - (x(i ,0,0) - x(i-1,0,0)) * sx(i-1,0,0); + return dxinv[0]*dxinv[0] * y; + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -Real mlndlap_adotx_aa (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&, Array4 const&, - GpuArray const&) noexcept -{ return Real(0.0); } +Real mlndlap_adotx_aa (int i, int j, int k, Array4 const& x, + Array4 const& sx, Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + return mlndlap_adotx_ha(i,j,k,x,sx,msk,dxinv); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_normalize_ha (int /*i*/, int /*j*/, int /*k*/, Array4 const&, Array4 const&, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_normalize_ha (int i, int, int, Array4 const& x, + Array4 const& sx, Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + if (!msk(i,0,0)) { + x(i,0,0) /= -dxinv[0]*dxinv[0] * (sx(i-1,0,0)+sx(i,0,0)); + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_normalize_aa (int /*i*/, int /*j*/, int /*k*/, Array4 const&, Array4 const&, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_normalize_aa (int i, int j, int k, Array4 const& x, + Array4 const& sx, Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + mlndlap_normalize_ha(i,j,k,x,sx,msk,dxinv); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_jacobi_ha (int, int, int, Array4 const&, Real, - Array4 const&, Array4 const&, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_jacobi_ha (int i, int, int, Array4 const& sol, Real Ax, + Array4 const& rhs, Array4 const& sx, + Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + if (msk(i,0,0)) { + sol(i,0,0) = Real(0.0); + } else { + sol(i,0,0) += Real(2.0/3.0) * (rhs(i,0,0) - Ax) + / (-dxinv[0]*dxinv[0] * (sx(i-1,0,0)+sx(i,0,0))); + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_jacobi_ha (Box const&, Array4 const&, Array4 const&, - Array4 const&, Array4 const&, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_jacobi_ha (Box const& bx, Array4 const& sol, Array4 const& Ax, + Array4 const& rhs, Array4 const& sx, + Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + Real fac = -dxinv[0]*dxinv[0]; + + amrex::LoopConcurrent(bx, [=] (int i, int, int) noexcept + { + if (msk(i,0,0)) { + sol(i,0,0) = Real(0.0); + } else { + sol(i,0,0) += Real(2.0/3.0) * (rhs(i,0,0) - Ax(i,0,0)) + / (fac * (sx(i-1,0,0)+sx(i,0,0))); + } + }); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_jacobi_aa (int, int, int, Array4 const&, Real, - Array4 const&, Array4 const&, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_jacobi_aa (int i, int j, int k, Array4 const& sol, Real Ax, + Array4 const& rhs, Array4 const& sig, + Array4 const& msk, GpuArray const& dxinv) noexcept +{ + mlndlap_jacobi_ha(i,j,k,sol,Ax,rhs,sig,msk,dxinv); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_jacobi_aa (Box const&, Array4 const&, Array4 const&, - Array4 const&, Array4 const&, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_jacobi_aa (Box const& bx, Array4 const& sol, Array4 const& Ax, + Array4 const& rhs, Array4 const& sig, + Array4 const& msk, GpuArray const& dxinv) noexcept +{ + mlndlap_jacobi_ha(bx,sol,Ax,rhs,sig,msk,dxinv); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_jacobi_c (int, int, int, Array4 const&, Real, - Array4 const&, Real, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_jacobi_c (int i, int, int, Array4 const& sol, Real Ax, + Array4 const& rhs, Real sig, + Array4 const& msk, GpuArray const& dxinv) noexcept +{ + if (msk(i,0,0)) { + sol(i,0,0) = Real(0.0); + } else { + sol(i,0,0) += Real(2.0/3.0) * (rhs(i,0,0) - Ax) + / (-dxinv[0]*dxinv[0]*Real(2.0)*sig); + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_jacobi_c (Box const&, Array4 const&, Array4 const&, - Array4 const&, Real, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_jacobi_c (Box const& bx, Array4 const& sol, Array4 const& Ax, + Array4 const& rhs, Real sig, + Array4 const& msk, GpuArray const& dxinv) noexcept +{ + amrex::LoopConcurrent(bx, [=] (int i, int, int) noexcept + { + if (msk(i,0,0)) { + sol(i,0,0) = Real(0.0); + } else { + sol(i,0,0) += Real(2.0/3.0) * (rhs(i,0,0) - Ax(i,0,0)) + / (-dxinv[0]*dxinv[0]*Real(2.0)*sig); + } + }); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_gauss_seidel_ha (Box const&, Array4 const&, - Array4 const&, Array4 const&, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_gauss_seidel_ha (Box const& bx, Array4 const& sol, + Array4 const& rhs, + Array4 const& sx, + Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + Real fac = dxinv[0]*dxinv[0]; + + amrex::Loop(bx, [=] (int i, int, int) noexcept + { + if (msk(i,0,0)) { + sol(i,0,0) = Real(0.0); + } else { + Real s0 = Real(-1.0) * fac * (sx(i-1,0,0)+sx(i,0,0)); + Real Ax = sol(i-1,0,0)*fac*sx(i-1,0,0) + + sol(i+1,0,0)*fac*sx(i ,0,0) + + sol(i ,0,0)*s0; + sol(i,0,0) += (rhs(i,0,0) - Ax) / s0; + } + }); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_gauss_seidel_aa (Box const&, Array4 const&, - Array4 const&, Array4 const&, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_gauss_seidel_aa (Box const& bx, Array4 const& sol, + Array4 const& rhs, + Array4 const& sx, + Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + mlndlap_gauss_seidel_ha(bx,sol,rhs,sx,msk,dxinv); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_gauss_seidel_c (Box const&, Array4 const&, - Array4 const&, Real, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_gauss_seidel_c (Box const& bx, Array4 const& sol, + Array4 const& rhs, Real sig, + Array4 const& msk, + GpuArray const& dxinv) noexcept +{ + Real fac = dxinv[0]*dxinv[0]; + + amrex::Loop(bx, [=] (int i, int, int) noexcept + { + if (msk(i,0,0)) { + sol(i,0,0) = Real(0.0); + } else { + Real s0 = Real(-2.0) * fac * sig; + Real Ax = sol(i-1,0,0)*fac*sig + + sol(i+1,0,0)*fac*sig + + sol(i ,0,0)*s0; + sol(i,0,0) += (rhs(i,0,0) - Ax) / s0; + } + }); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_gauss_seidel_with_line_solve_aa(Box const&, Array4 const&, Array4 const&, Array4 const&, Array4 const&, GpuArray const&) noexcept -{} +{ + amrex::Abort("mlndlap_gauss_seidel_with_line_solve_aa: not implemented in 1D"); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_restriction (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&, Array4 const&) noexcept -{} +void mlndlap_restriction (int i, int, int, Array4 const& crse, + Array4 const& fine, Array4 const& msk) noexcept +{ + int ii = i*2; + if (msk(ii,0,0)) { + crse(i,0,0) = Real(0.0); + } else { + crse(i,0,0) = Real(1./4.) *(fine(ii-1,0,0) + + Real(2.)* fine(ii ,0,0) + + fine(ii+1,0,0)); + } +} template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_restriction (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&, Array4 const&, - Box const&, - GpuArray const&, - GpuArray const&) noexcept -{} +void mlndlap_restriction (int i, int, int, Array4 const& crse, + Array4 const& fine, Array4 const& msk, + Box const& fdom, + GpuArray const& bclo, + GpuArray const& bchi) noexcept + +{ + const int ii = i*rr; + if (msk(ii,0,0)) { + crse(i,0,0) = Real(0.0); + } else { + const auto ndlo = fdom.smallEnd(0); + const auto ndhi = fdom.bigEnd(0); + Real tmp = Real(0.0); + for (int ioff = -rr+1; ioff <= rr-1; ++ioff) { + Real wx = rr - std::abs(ioff); + int itmp = ii + ioff; + if ((itmp < ndlo && (bclo[0] == LinOpBCType::Neumann || + bclo[0] == LinOpBCType::inflow)) || + (itmp > ndhi && (bchi[0] == LinOpBCType::Neumann || + bchi[0] == LinOpBCType::inflow))) { + itmp = ii - ioff; + } + tmp += wx*fine(itmp,0,0); + } + crse(i,0,0) = tmp*(Real(1.0)/Real(rr*rr)); + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_semi_restriction (int /*i*/, int /*j*/, int /*k*/, Array4 const&, Array4 const&, Array4 const&, int) noexcept -{} +{ + amrex::Abort("mlndlap_semi_restriction: not implemented in 1D"); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_interpadd_c (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&, - Array4 const&) noexcept -{} +void mlndlap_interpadd_c (int i, int , int, Array4 const& fine, + Array4 const& crse, + Array4 const& msk) noexcept +{ + if (!msk(i,0,0)) { + int ic = amrex::coarsen(i,2); + bool i_is_odd = (ic*2 != i); + if (i_is_odd) { + // Node on X line + fine(i,0,0) += Real(0.5) * (crse(ic,0,0)+crse(ic+1,0,0)); + } else { + // Node coincident with coarse node + fine(i,0,0) += crse(ic,0,0); + } + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_interpadd_aa (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&, Array4 const&, - Array4 const&) noexcept -{} +void mlndlap_interpadd_aa (int i, int , int, Array4 const& fine, + Array4 const& crse, Array4 const& sig, + Array4 const& msk) noexcept +{ + if (!msk(i,0,0)) { + int ic = amrex::coarsen(i,2); + bool i_is_odd = (ic*2 != i); + if (i_is_odd) { + // Node on X line + fine(i,0,0) += (sig(i-1,0,0)*crse(ic,0,0) + sig(i,0,0)*crse(ic+1,0,0)) + / (sig(i-1,0,0) + sig(i,0,0)); + } else { + // Node coincident with coarse node + fine(i,0,0) += crse(ic,0,0); + } + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_semi_interpadd_aa (int /*i*/, int /*j*/, int /*k*/, Array4 const&, @@ -171,34 +432,59 @@ void mlndlap_semi_interpadd_aa (int /*i*/, int /*j*/, int /*k*/, Array4 co {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_interpadd_ha (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&, Array4 const&, - Array4 const&) noexcept -{} +void mlndlap_interpadd_ha (int i, int j, int k, Array4 const& fine, + Array4 const& crse, Array4 const& sig, + Array4 const& msk) noexcept +{ + mlndlap_interpadd_aa(i,j,k,fine,crse,sig,msk); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_divu (int /*i*/, int /*j*/, int /*k*/, Array4 const&, Array4 const&, - Array4 const&, - GpuArray const&, +void mlndlap_divu (int i, int, int, Array4 const& rhs, Array4 const& vel, + Array4 const& msk, GpuArray const& dxinv, Box const&, GpuArray const&, GpuArray const&) noexcept -{} +{ + Real fac = dxinv[0]; + if (msk(i,0,0)) { + rhs(i,0,0) = Real(0.0); + } else { + rhs(i,0,0) = fac*(vel(i,0,0)-vel(i-1,0,0)); + } +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -Real mlndlap_rhcc (int /*i*/, int /*j*/, int /*k*/, Array4 const&, - Array4 const&) noexcept -{ return Real(0.); } +Real mlndlap_rhcc (int i, int, int, Array4 const& rhcc, + Array4 const& msk) noexcept +{ + Real r; + if (msk(i,0,0)) { + r = Real(0.0); + } else { + r = Real(0.5) * (rhcc(i-1,0,0)+rhcc(i,0,0)); + } + return r; +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_mknewu (int /*i*/, int /*j*/, int /*k*/, Array4 const&, Array4 const&, - Array4 const&, GpuArray const&) noexcept -{} +void mlndlap_mknewu (int i, int, int, Array4 const& u, + Array4 const& p, + Array4 const& sig, + GpuArray const& dxinv) noexcept +{ + Real fac = dxinv[0]; + u(i,0,0) -= sig(i,0,0)*fac*(p(i+1,0,0)-p(i,0,0)); +} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void mlndlap_mknewu_c (int /*i*/, int /*j*/, int /*k*/, Array4 const&, Array4 const&, - Real, GpuArray const&) noexcept -{} +void mlndlap_mknewu_c (int i, int, int, Array4 const& u, + Array4 const& p, + Real sig, GpuArray const& dxinv) noexcept +{ + Real fac = dxinv[0]; + u(i,0,0) -= sig*fac*(p(i+1,0,0)-p(i,0,0)); +} template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_misc.cpp b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_misc.cpp index 9cb5ec880f..2c3e7bc342 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_misc.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_misc.cpp @@ -113,14 +113,19 @@ MLNodeLaplacian::averageDownCoeffsSameAmrLevel (int amrlev) for (int mglev = 1; mglev < m_num_mg_levels[amrlev]; ++mglev) { + IntVect ratio = mg_coarsen_ratio_vec[mglev-1]; +#if (AMREX_SPACEDIM == 1) + int idir = 0; + bool regular_coarsening = true; +#else int idir = 2; bool regular_coarsening = mg_coarsen_ratio_vec[mglev-1] == mg_coarsen_ratio; - IntVect ratio = mg_coarsen_ratio_vec[mglev-1]; if (ratio[1] == 1) { idir = 1; } else if (ratio[0] == 1) { idir = 0; } +#endif for (int idim = 0; idim < nsigma; ++idim) { const MultiFab& fine = *m_sigma[amrlev][mglev-1][idim]; @@ -891,6 +896,7 @@ MLNodeLaplacian::compRHS (const Vector& rhs, const Vector& { #if (AMREX_SPACEDIM == 1) amrex::ignore_unused(rhs,vel,rhnd,a_rhcc); + amrex::Abort("MLNodeLaplacian::compRHS: 1D not supported"); #else // // Note that div vel we copmute on a coarse/fine nodes is not a diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_sync.cpp b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_sync.cpp index 587f9508ef..07cf1628f4 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_sync.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_sync.cpp @@ -19,6 +19,10 @@ MLNodeLaplacian::compSyncResidualCoarse (MultiFab& sync_resid, const MultiFab& a { BL_PROFILE("MLNodeLaplacian::SyncResCrse()"); +#if (AMREX_SPACEDIM == 1) + amrex::Abort("MLNodeLaplacian::compSyncResidualCoarse: 1D not supported"); +#endif + sync_resid.setVal(0.0); const Geometry& geom = m_geom[0][0]; @@ -337,6 +341,10 @@ MLNodeLaplacian::compSyncResidualFine (MultiFab& sync_resid, const MultiFab& phi { BL_PROFILE("MLNodeLaplacian::SyncResFine()"); +#if (AMREX_SPACEDIM == 1) + amrex::Abort("MLNodeLaplacian::compSyncResidualFine: 1D not supported"); +#endif + const auto& sigma_orig = m_sigma[0][0][0]; const iMultiFab& dmsk = *m_dirichlet_mask[0][0]; @@ -617,6 +625,7 @@ MLNodeLaplacian::reflux (int crse_amrlev, { #if (AMREX_SPACEDIM == 1) amrex::ignore_unused(crse_amrlev,res,crse_sol,crse_rhs,a_fine_res,fine_sol,fine_rhs); + amrex::Abort("MLNodeLaplacian::reflux: 1D not supported"); #else // // Note that the residue we copmute on a coarse/fine node is not a diff --git a/Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp b/Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp index 9a48922a60..9999fe7d12 100644 --- a/Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp +++ b/Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp @@ -154,18 +154,16 @@ MyTest::solveABecLaplacian () mlabec.setMaxOrder(linop_maxorder); - // This is a 3d problem with homogeneous Neumann BC - mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann, + mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Dirichlet, LinOpBCType::Neumann, LinOpBCType::Neumann)}, {AMREX_D_DECL(LinOpBCType::Neumann, - LinOpBCType::Neumann, + LinOpBCType::Dirichlet, LinOpBCType::Neumann)}); for (int ilev = 0; ilev < nlevels; ++ilev) { - // for problem with pure homogeneous Neumann BC, we could pass a nullptr - mlabec.setLevelBC(ilev, nullptr); + mlabec.setLevelBC(ilev, &solution[ilev]); } mlabec.setScalars(ascalar, bscalar); @@ -213,20 +211,18 @@ MyTest::solveABecLaplacian () mlabec.setMaxOrder(linop_maxorder); - // This is a 3d problem with homogeneous Neumann BC - mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann, + mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Dirichlet, LinOpBCType::Neumann, LinOpBCType::Neumann)}, {AMREX_D_DECL(LinOpBCType::Neumann, - LinOpBCType::Neumann, + LinOpBCType::Dirichlet, LinOpBCType::Neumann)}); if (ilev > 0) { mlabec.setCoarseFineBC(&solution[ilev-1], ref_ratio); } - // for problem with pure homogeneous Neumann BC, we could pass a nullptr - mlabec.setLevelBC(0, nullptr); + mlabec.setLevelBC(0, &solution[ilev]); mlabec.setScalars(ascalar, bscalar); @@ -263,18 +259,6 @@ MyTest::solveABecLaplacian () mlmg.solve({&solution[ilev]}, {&rhs[ilev]}, tol_rel, tol_abs); } } - - // Since this problem has Neumann BC, solution + constant is also a - // solution. So we are going to shift the solution by a constant - // for comparison with the "exact solution". - // The statement above is incorrect because we have the a term, albeit small. - const Real npts = grids[0].d_numPts(); - const Real avg1 = exact_solution[0].sum(); - const Real avg2 = solution[0].sum(); - const Real offset = (avg1-avg2)/npts; - for (int ilev = 0; ilev < nlevels; ++ilev) { - solution[ilev].plus(offset, 0, 1, 0); - } } void @@ -406,18 +390,6 @@ MyTest::solveABecLaplacianInhomNeumann () mlmg.solve({&solution[ilev]}, {&rhs[ilev]}, tol_rel, tol_abs); } } - - // Since this problem has Neumann BC, solution + constant is also a - // solution. So we are going to shift the solution by a constant - // for comparison with the "exact solution". - // The statement above is incorrect because we have the a term, albeit small. - const Real npts = grids[0].d_numPts(); - const Real avg1 = exact_solution[0].sum(); - const Real avg2 = solution[0].sum(); - const Real offset = (avg1-avg2)/npts; - for (int ilev = 0; ilev < nlevels; ++ilev) { - solution[ilev].plus(offset, 0, 1, 0); - } } void @@ -481,23 +453,22 @@ MyTest::solveABecLaplacianGMRES () const auto tol_rel = Real(1.e-10); const auto tol_abs = Real(0.0); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(composite_solve == false, - "solveABecLaplacianGMRES does not support composite solve"); - const auto nlevels = static_cast(geom.size()); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(composite_solve == false || nlevels == 1, + "solveABecLaplacianGMRES does not support composite solve"); + for (int ilev = 0; ilev < nlevels; ++ilev) { MLABecLaplacian mlabec({geom[ilev]}, {grids[ilev]}, {dmap[ilev]}, info); mlabec.setMaxOrder(linop_maxorder); - // This is a 3d problem with homogeneous Neumann BC - mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Neumann, + mlabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Dirichlet, LinOpBCType::Neumann, LinOpBCType::Neumann)}, {AMREX_D_DECL(LinOpBCType::Neumann, - LinOpBCType::Neumann, + LinOpBCType::Dirichlet, LinOpBCType::Neumann)}); if (ilev > 0) { @@ -505,7 +476,7 @@ MyTest::solveABecLaplacianGMRES () } // for problem with pure homogeneous Neumann BC, we could pass a nullptr - mlabec.setLevelBC(0, nullptr); + mlabec.setLevelBC(0, &solution[ilev]); mlabec.setScalars(ascalar, bscalar); @@ -536,18 +507,6 @@ MyTest::solveABecLaplacianGMRES () << " " << res.norm1(0) << " " << res.norm2(0) << '\n'; } } - - // Since this problem has Neumann BC, solution + constant is also a - // solution. So we are going to shift the solution by a constant - // for comparison with the "exact solution". - // The statement above is incorrect because we have the a term, albeit small. - const Real npts = grids[0].d_numPts(); - const Real avg1 = exact_solution[0].sum(); - const Real avg2 = solution[0].sum(); - const Real offset = (avg1-avg2)/npts; - for (int ilev = 0; ilev < nlevels; ++ilev) { - solution[ilev].plus(offset, 0, 1, 0); - } } void diff --git a/Tests/LinearSolvers/ABecLaplacian_C/initProb.cpp b/Tests/LinearSolvers/ABecLaplacian_C/initProb.cpp index ec102ae060..59a9dbcdd2 100644 --- a/Tests/LinearSolvers/ABecLaplacian_C/initProb.cpp +++ b/Tests/LinearSolvers/ABecLaplacian_C/initProb.cpp @@ -45,10 +45,11 @@ MyTest::initProbABecLaplacian () #endif for (MFIter mfi(rhs[ilev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.tilebox(); + const Box& vbx = mfi.validbox(); const Box& gbx = mfi.growntilebox(1); auto rhsfab = rhs[ilev].array(mfi); + auto solfab = solution[ilev].array(mfi); auto exactfab = exact_solution[ilev].array(mfi); auto acoeffab = acoef[ilev].array(mfi); auto bcoeffab = bcoef[ilev].array(mfi); @@ -56,18 +57,10 @@ MyTest::initProbABecLaplacian () amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - actual_init_bcoef(i,j,k,bcoeffab,prob_lo,prob_hi,dx); - }); - - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - actual_init_abeclap(i,j,k,rhsfab,exactfab,acoeffab,bcoeffab, - a,b,prob_lo,prob_hi,dx); + actual_init_abeclap(i,j,k,rhsfab,solfab,exactfab,acoeffab,bcoeffab, + a,b,prob_lo,prob_hi,dx,vbx); }); } - - solution[ilev].setVal(0.0); } } diff --git a/Tests/LinearSolvers/ABecLaplacian_C/initProb_K.H b/Tests/LinearSolvers/ABecLaplacian_C/initProb_K.H index 0b6fbe20f7..8d98aa1078 100644 --- a/Tests/LinearSolvers/ABecLaplacian_C/initProb_K.H +++ b/Tests/LinearSolvers/ABecLaplacian_C/initProb_K.H @@ -2,6 +2,7 @@ #define INIT_PROB_K_H_ #include +#include AMREX_GPU_DEVICE AMREX_FORCE_INLINE void actual_init_poisson (int i, int j, int k, @@ -70,14 +71,16 @@ void actual_init_bcoef (int i, int j, int k, AMREX_GPU_DEVICE AMREX_FORCE_INLINE void actual_init_abeclap (int i, int j, int k, - amrex::Array4 const& rhs, - amrex::Array4 const& exact, - amrex::Array4 const& alpha, - amrex::Array4 const& beta, + amrex::Array4 const& rhs, + amrex::Array4 const& sol, + amrex::Array4 const& exact, + amrex::Array4 const& alpha, + amrex::Array4 const& beta, amrex::Real a, amrex::Real b, amrex::GpuArray const& prob_lo, amrex::GpuArray const& prob_hi, - amrex::GpuArray const& dx) + amrex::GpuArray const& dx, + amrex::Box const& vbx) { constexpr amrex::Real w = 0.05; constexpr amrex::Real sigma = 10.; @@ -105,25 +108,43 @@ void actual_init_abeclap (int i, int j, int k, #endif amrex::Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); - amrex::Real tmp = std::cosh(theta*(r-0.25)); - amrex::Real dbdrfac = (sigma-1.)/2./(tmp*tmp) * theta/r; - dbdrfac *= b; - - alpha(i,j,k) = 1.; - exact(i,j,k) = std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z) - + .25 * std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z); + beta(i,j,k) = (sigma-1.)/2.*std::tanh(theta*(r-0.25)) + (sigma+1.)/2.; - rhs(i,j,k) = beta(i,j,k)*b*fac*(std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z) - + std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z)) - + dbdrfac*((x-xc)*(tpi*std::sin(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z) - + pi*std::sin(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z)) - + (y-yc)*(tpi*std::cos(tpi*x) * std::sin(tpi*y) * std::cos(tpi*z) - + pi*std::cos(fpi*x) * std::sin(fpi*y) * std::cos(fpi*z)) - + (z-zc)*(tpi*std::cos(tpi*x) * std::cos(tpi*y) * std::sin(tpi*z) - + pi*std::cos(fpi*x) * std::cos(fpi*y) * std::sin(fpi*z))) - + a * (std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z) - + 0.25 * std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z)); + if (vbx.contains(i,j,k)) { + amrex::Real tmp = std::cosh(theta*(r-0.25)); + amrex::Real dbdrfac = (sigma-1.)/2./(tmp*tmp) * theta/r; + dbdrfac *= b; + + alpha(i,j,k) = 1.; + + sol(i,j,k) = 0.; // Use zero as initial guess for the linear solver + + exact(i,j,k) = std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z) + + .25 * std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z); + + rhs(i,j,k) = beta(i,j,k)*b*fac*(std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z) + + std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z)) + + dbdrfac*((x-xc)*(tpi*std::sin(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z) + + pi*std::sin(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z)) + + (y-yc)*(tpi*std::cos(tpi*x) * std::sin(tpi*y) * std::cos(tpi*z) + + pi*std::cos(fpi*x) * std::sin(fpi*y) * std::cos(fpi*z)) + + (z-zc)*(tpi*std::cos(tpi*x) * std::cos(tpi*y) * std::sin(tpi*z) + + pi*std::cos(fpi*x) * std::cos(fpi*y) * std::sin(fpi*z))) + + a * (std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z) + + 0.25 * std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z)); + } else { + amrex::Real xb = std::clamp(x, prob_lo[0], prob_hi[0]); + amrex::Real yb = std::clamp(y, prob_lo[1], prob_hi[1]); +#if (AMREX_SPACEDIM == 2) + amrex::Real zb = 0.0; +#else + amrex::Real zb = std::clamp(z, prob_lo[2], prob_hi[2]); +#endif + // provide Dirichlet BC + sol(i,j,k) = std::cos(tpi*xb) * std::cos(tpi*yb) * std::cos(tpi*zb) + + .25 * std::cos(fpi*xb) * std::cos(fpi*yb) * std::cos(fpi*zb); + } } AMREX_GPU_DEVICE AMREX_FORCE_INLINE diff --git a/Tests/LinearSolvers/NodalPoisson/MyTest.cpp b/Tests/LinearSolvers/NodalPoisson/MyTest.cpp index ebb73dc123..39c781d879 100644 --- a/Tests/LinearSolvers/NodalPoisson/MyTest.cpp +++ b/Tests/LinearSolvers/NodalPoisson/MyTest.cpp @@ -251,8 +251,16 @@ MyTest::initData () constexpr Real fac = tpi*tpi*AMREX_SPACEDIM; Real x = i*dx[0]; +#if (AMREX_SPACEDIM > 1) Real y = j*dx[1]; +#else + Real y = Real(0.0); +#endif +#if (AMREX_SPACEDIM > 2) Real z = k*dx[2]; +#else + Real z = Real(0.0); +#endif phi(i,j,k) = (std::cos(tpi*x) * std::cos(tpi*y) * std::cos(tpi*z)) + 0.25 * (std::cos(fpi*x) * std::cos(fpi*y) * std::cos(fpi*z));