diff --git a/Docs/sphinx_documentation/source/LinearSolvers.rst b/Docs/sphinx_documentation/source/LinearSolvers.rst index 21fa7fb95e..d4a3e971f4 100644 --- a/Docs/sphinx_documentation/source/LinearSolvers.rst +++ b/Docs/sphinx_documentation/source/LinearSolvers.rst @@ -863,14 +863,16 @@ discretized form of .. math:: \nabla \times (\alpha \nabla \times \vec{E}) + \beta \vec{E} = \vec{f}, -where :math:`\vec{E}` and :math:`\vec{f}` are defined on cell edges, -:math:`\alpha` is a positive scalar, and :math:`\beta` is a non-negative -scalar (either a constant or a field). An :cpp:`Array` of three -:cpp:`MultiFab`\ s is used to store the components of :math:`\vec{E}` and -:math:`\vec{f}`. It is the user's responsibility to ensure that the -right-hand-side data are consistent on edges shared by multiple -:cpp:`Box`\ es. If needed, you can call :cpp:`MLCurlCurl::prepareRHS` to -perform this synchronization. +where :math:`\vec{E}` and :math:`\vec{f}` are defined on cell edges. The +coefficient :math:`\alpha` may be supplied either as a single positive +scalar through :cpp:`MLCurlCurl::setScalars`, or as a nodal :cpp:`MultiFab` +by calling :cpp:`MLCurlCurl::setAlpha` with one entry per AMR level. The +:math:`\beta` term can be a non-negative scalar or an edge-centered field +set via :cpp:`setBeta`. An :cpp:`Array` of three :cpp:`MultiFab`\ s is used +to store the components of :math:`\vec{E}` and :math:`\vec{f}`. It is the +user's responsibility to ensure that the right-hand-side data are consistent +on edges shared by multiple :cpp:`Box`\ es. If needed, you can call +:cpp:`MLCurlCurl::prepareRHS` to perform this synchronization. The solver supports 1D, 2D and 3D. Note that even in the 1D and 2D cases, :math:`\vec{E}` still has three components, one for each spatial diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H index 456c135d05..78e7a8bf27 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H @@ -10,8 +10,10 @@ 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 either a non-negative scalar or a MultiFab. + * Here E is an Array of 3 MultiFabs on the staggered grid. The coefficient + * \f$\alpha\f$ can be supplied either as a positive scalar (via setScalars) + * or as a nodal MultiFab (via setAlpha), and \f$\beta\f$ can be either a + * non-negative scalar or an edge-centered MultiFab (via setBeta). * * It's the caller's responsibility to make sure rhs has consistent nodal * data. If needed, one could call prepareRHS for this. @@ -45,11 +47,21 @@ public: const LPInfo& a_info = LPInfo(), int a_coord = 0); + //! Set scalar coefficients. This clears any previously supplied nodal alpha + //! or edge-centered beta MultiFabs, so call setAlpha/setBeta afterwards if + //! you still need spatially varying coefficients. void setScalars (RT a_alpha, RT a_beta) noexcept; - //! This is needed only if there is variable beta coefficient. + //! This is needed only if there is variable beta coefficient. Call this + //! after setScalars if both APIs are used, because setScalars clears the + //! cached beta MultiFabs. void setBeta (const Vector>& a_bcoefs); + //! Set variable alpha defined on nodal MultiFab. Call this after + //! setScalars if both APIs are used, because setScalars clears the cached + //! nodal alpha fields. + void setAlpha (const Vector& a_acoeffs); + //! Synchronize RHS on nodal points. If the user can guarantee it, this //! function does not need to be called. void prepareRHS (Vector const& rhs) const; @@ -140,6 +152,8 @@ private: void update_lusolver (); void set_curvilinear_domain_bc (); + void clearAlphaMultiFab (); + void clearBetaMultiFab (); int m_coord = 0; @@ -160,6 +174,7 @@ private: Vector>>>> m_lusolver; Vector,3>>> m_bcoefs; + Vector,3>>> m_acoefs; bool m_use_pcg = false; bool m_needs_update = true; }; diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp index 59c7c33e52..44257ccf0a 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp @@ -1,5 +1,7 @@ #include #include +#include +#include namespace amrex { @@ -40,6 +42,11 @@ void MLCurlCurl::define (const Vector& a_geom, m_bcoefs[amrlev].resize(this->m_num_mg_levels[amrlev]); } + m_acoefs.resize(this->m_num_amr_levels); + for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) { + m_acoefs[amrlev].resize(this->m_num_mg_levels[amrlev]); + } + m_lusolver.resize(this->m_num_amr_levels); for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) { m_lusolver[amrlev].resize(this->m_num_mg_levels[amrlev]); @@ -51,6 +58,8 @@ void MLCurlCurl::setScalars (RT a_alpha, RT a_beta) noexcept m_needs_update = true; m_alpha = a_alpha; m_beta = a_beta; + clearAlphaMultiFab(); + clearBetaMultiFab(); AMREX_ASSERT(m_beta > RT(0)); } @@ -114,6 +123,108 @@ void MLCurlCurl::setBeta (const Vector>& a_bcoefs) #endif } } + + for (auto& amrvec : m_lusolver) { + for (auto& mgptr : amrvec) { + mgptr.reset(); + } + } +} + +void MLCurlCurl::setAlpha (const Vector& a_acoeffs) +{ + AMREX_ALWAYS_ASSERT(static_cast(a_acoeffs.size()) == m_num_amr_levels); + + m_needs_update = true; + + auto lobc = LoBC(); + auto hibc = HiBC(); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (lobc[idim] != LinOpBCType::Periodic) { lobc[idim] = LinOpBCType::Neumann; } + if (hibc[idim] != LinOpBCType::Periodic) { hibc[idim] = LinOpBCType::Neumann; } + } + + for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) { + AMREX_ALWAYS_ASSERT(a_acoeffs[amrlev]->is_nodal()); + MultiFab nodal_alpha(a_acoeffs[amrlev]->boxArray(), + a_acoeffs[amrlev]->DistributionMap(), 1, 1); + MultiFab::Copy(nodal_alpha, *a_acoeffs[amrlev], 0, 0, 1, 0); + nodal_alpha.FillBoundaryAndSync(m_geom[amrlev][0].periodicity()); + + for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev) { + + Box nd_domain = amrex::surroundingNodes(m_geom[amrlev][mglev].Domain()); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(nodal_alpha); mfi.isValid(); ++mfi) { + auto const& afab = nodal_alpha.array(mfi); + Box const& box = mfi.validbox(); + mlndlap_applybc(box, afab, nd_domain, lobc, hibc); + } + + auto const& anode = nodal_alpha.const_arrays(); + + GpuArray,3> aface; + for (int idim = 0; idim < 3; ++idim) { + IntVect typ(0); + if (idim < AMREX_SPACEDIM) { typ[idim] = 1; } + m_acoefs[amrlev][mglev][idim] = std::make_unique + (amrex::convert(m_grids[amrlev][mglev], typ), + m_dmap[amrlev][mglev], 1, 1); + aface[idim] = m_acoefs[amrlev][mglev][idim]->arrays(); + } + + amrex::ParallelFor(nodal_alpha, IntVect(1), + [=] AMREX_GPU_DEVICE (int b, int i, int j, int k) + { + auto const& an = anode[b]; + auto const& ax = aface[0][b]; + auto const& ay = aface[1][b]; + auto const& az = aface[2][b]; + if (ax.contains(i,j,k)) { +#if (AMREX_SPACEDIM == 1) + ax(i,0,0) = an(i,0,0); +#elif (AMREX_SPACEDIM == 2) + ax(i,j,0) = Real(0.5)*(an(i,j,0)+an(i,j+1,0)); +#else + ax(i,j,k) = Real(0.25)*(an(i,j,k)+an(i,j+1,k)+an(i,j,k+1)+an(i,j+1,k+1)); +#endif + } + if (ay.contains(i,j,k)) { +#if (AMREX_SPACEDIM == 1) + ay(i,0,0) = Real(0.5)*(an(i,0,0)+an(i+1,0,0)); +#elif (AMREX_SPACEDIM == 2) + ay(i,j,0) = Real(0.5)*(an(i,j,0)+an(i+1,j,0)); +#else + ay(i,j,k) = Real(0.25)*(an(i,j,k)+an(i+1,j,k)+an(i,j,k+1)+an(i+1,j,k+1)); +#endif + } + if (az.contains(i,j,k)) { +#if (AMREX_SPACEDIM == 1) + az(i,j,k) = an(i,j,k); +#elif (AMREX_SPACEDIM >= 2) + az(i,j,k) = Real(0.25)*(an(i,j,k)+an(i+1,j,k)+an(i,j+1,k)+an(i+1,j+1,k)); +#endif + } + }); + + if (mglev+1 < m_num_mg_levels[amrlev]) { + MultiFab tmp(amrex::convert(m_grids[amrlev][mglev+1], IntVect(1)), + m_dmap[amrlev][mglev+1], 1, 1); + IntVect ratio = (amrlev > 0) ? IntVect(2) : mg_coarsen_ratio_vec[mglev]; + average_down_nodal(nodal_alpha, tmp, ratio); + std::swap(nodal_alpha, tmp); + nodal_alpha.FillBoundary(m_geom[amrlev][mglev+1].periodicity()); + } + } + } + + for (auto& amrvec : m_lusolver) { + for (auto& mgptr : amrvec) { + mgptr.reset(); + } + } } void MLCurlCurl::prepareRHS (Vector const& rhs) const @@ -265,11 +376,14 @@ MLCurlCurl::apply (int amrlev, int mglev, MF& out, MF& in, BCMode /*bc_mode*/, { applyBC(amrlev, mglev, in, CurlCurlStateType::x); - auto adxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto adxinv = dxinv; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { adxinv[idim] *= std::sqrt(m_alpha); } auto const b = m_beta; + bool const has_beta = (m_bcoefs[amrlev][mglev][0] != nullptr); + bool const has_alpha = (m_acoefs[amrlev][mglev][0] != nullptr); auto dinfo = getDirichletInfo(amrlev,mglev); auto coord = m_coord; @@ -289,7 +403,57 @@ MLCurlCurl::apply (int amrlev, int mglev, MF& out, MF& in, BCMode /*bc_mode*/, auto const& xin = in[0].array(mfi); auto const& yin = in[1].array(mfi); auto const& zin = in[2].array(mfi); - if (m_bcoefs[amrlev][mglev][0]) { + + if (has_alpha) { + Array4 bcx, bcy, bcz; + if (has_beta) { + bcx = m_bcoefs[amrlev][mglev][0]->const_array(mfi); + bcy = m_bcoefs[amrlev][mglev][1]->const_array(mfi); + bcz = m_bcoefs[amrlev][mglev][2]->const_array(mfi); + } + auto const afx = m_acoefs[amrlev][mglev][0]->const_array(mfi); + auto const afy = m_acoefs[amrlev][mglev][1]->const_array(mfi); + auto const afz = m_acoefs[amrlev][mglev][2]->const_array(mfi); + amrex::ParallelFor(xbx, ybx, zbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (dinfo.is_dirichlet_x_edge(i,j,k)) { + xout(i,j,k) = Real(0.0); + } else { + Real beta = bcx ? bcx(i,j,k) : b; + mlcurlcurl_adotx_x_alpha(i,j,k,xout,xin,yin,zin, + afy,afz,beta,dxinv); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (dinfo.is_dirichlet_y_edge(i,j,k)) { + yout(i,j,k) = Real(0.0); + } else { + Real beta = bcy ? bcy(i,j,k) : b; + mlcurlcurl_adotx_y_alpha(i,j,k,yout,xin,yin,zin, + afx,afz,beta,dxinv +#if (AMREX_SPACEDIM < 3) + ,coord +#endif + ); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (dinfo.is_dirichlet_z_edge(i,j,k)) { + zout(i,j,k) = Real(0.0); + } else { + Real beta = bcz ? bcz(i,j,k) : b; + mlcurlcurl_adotx_z_alpha(i,j,k,zout,xin,yin,zin, + afx,afy,beta,dxinv +#if (AMREX_SPACEDIM < 3) + ,coord +#endif + ); + } + }); + } else if (has_beta) { auto const& bcx = m_bcoefs[amrlev][mglev][0]->const_array(mfi); auto const& bcy = m_bcoefs[amrlev][mglev][1]->const_array(mfi); auto const& bcz = m_bcoefs[amrlev][mglev][2]->const_array(mfi); @@ -405,7 +569,8 @@ void MLCurlCurl::smooth1D (int amrlev, int mglev, MF& sol, MF const& rhs, auto b = m_beta; auto dinfo = getDirichletInfo(amrlev,mglev); - auto adxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto adxinv = dxinv; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { adxinv[idim] *= std::sqrt(m_alpha); } @@ -417,7 +582,37 @@ void MLCurlCurl::smooth1D (int amrlev, int mglev, MF& sol, MF const& rhs, MultiFab nmf(amrex::convert(rhs[0].boxArray(),IntVect(1)), rhs[0].DistributionMap(), 1, 0, MFInfo().SetAlloc(false)); - if (m_bcoefs[amrlev][mglev][0]) { + bool const has_beta = (m_bcoefs[amrlev][mglev][0] != nullptr); + bool const has_alpha = (m_acoefs[amrlev][mglev][0] != nullptr); + + if (has_alpha && has_beta) { + auto const& acy = m_acoefs[amrlev][mglev][1]->const_arrays(); + auto const& acz = m_acoefs[amrlev][mglev][2]->const_arrays(); + auto const& bcx = m_bcoefs[amrlev][mglev][0]->const_arrays(); + auto const& bcy = m_bcoefs[amrlev][mglev][1]->const_arrays(); + auto const& bcz = m_bcoefs[amrlev][mglev][2]->const_arrays(); + ParallelFor( nmf, [=] AMREX_GPU_DEVICE(int bno, int i, int j, int k) + { + bool valid_x = i <= xhi; // x is cell-centered, not nodal + mlcurlcurl_smooth_1d_alpha_beta(i,j,k,ex[bno],ey[bno],ez[bno], + rhsx[bno],rhsy[bno],rhsz[bno], + bcx[bno],bcy[bno],bcz[bno], + dxinv,color,dinfo,valid_x,coord, + acy[bno],acz[bno]); + }); + } else if (has_alpha && !has_beta) { + auto const& acy = m_acoefs[amrlev][mglev][1]->const_arrays(); + auto const& acz = m_acoefs[amrlev][mglev][2]->const_arrays(); + ParallelFor( nmf, [=] AMREX_GPU_DEVICE(int bno, int i, int j, int k) + { + bool valid_x = i <= xhi; // x is cell-centered, not nodal + mlcurlcurl_smooth_1d_alpha(i,j,k,ex[bno],ey[bno],ez[bno], + rhsx[bno],rhsy[bno],rhsz[bno], + b, + dxinv,color,dinfo,valid_x,coord, + acy[bno],acz[bno]); + }); + } else if (!has_alpha && has_beta) { auto const& bcx = m_bcoefs[amrlev][mglev][0]->const_arrays(); auto const& bcy = m_bcoefs[amrlev][mglev][1]->const_arrays(); auto const& bcz = m_bcoefs[amrlev][mglev][2]->const_arrays(); @@ -455,21 +650,26 @@ void MLCurlCurl::smooth4 (int amrlev, int mglev, MF& sol, MF const& rhs, auto const& rhsy = rhs[1].const_arrays(); auto const& rhsz = rhs[2].const_arrays(); -#if (AMREX_SPACEDIM == 2) - auto b = m_beta; -#endif - - auto adxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); + auto adxinv = dxinv; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { adxinv[idim] *= std::sqrt(m_alpha); } + bool const has_beta = (m_bcoefs[amrlev][mglev][0] != nullptr); + bool const has_alpha = (m_acoefs[amrlev][mglev][0] != nullptr); + bool const use_pcg = m_use_pcg || has_alpha; + // We support LU solver with variable beta and scalar alpha. + auto dinfo = getDirichletInfo(amrlev,mglev); auto sinfo = getSymmetryInfo(amrlev,mglev); MultiFab nmf(amrex::convert(rhs[0].boxArray(),IntVect(1)), rhs[0].DistributionMap(), 1, 0, MFInfo().SetAlloc(false)); - if (m_lusolver[amrlev][mglev]) { + if (m_lusolver[amrlev][mglev] && !has_alpha && !has_beta) { +#if (AMREX_SPACEDIM == 2) + auto b = m_beta; +#endif auto* plusolver = m_lusolver[amrlev][mglev]->dataPtr(); ParallelFor(nmf, [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k) { @@ -480,14 +680,48 @@ void MLCurlCurl::smooth4 (int amrlev, int mglev, MF& sol, MF const& rhs, #endif adxinv,color,*plusolver,dinfo,sinfo); }); + } else if (has_alpha && has_beta) { + auto const& acx = m_acoefs[amrlev][mglev][0]->const_arrays(); + auto const& acy = m_acoefs[amrlev][mglev][1]->const_arrays(); + auto const& acz = m_acoefs[amrlev][mglev][2]->const_arrays(); + auto b = m_beta; + auto const& bcx = m_bcoefs[amrlev][mglev][0]->const_arrays(); + auto const& bcy = m_bcoefs[amrlev][mglev][1]->const_arrays(); + auto const& bcz = m_bcoefs[amrlev][mglev][2]->const_arrays(); + ParallelFor(nmf, [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k) + { + mlcurlcurl_gs4_alpha(i,j,k,ex[bno],ey[bno],ez[bno], + rhsx[bno],rhsy[bno],rhsz[bno], + dxinv,color, + acx[bno],acy[bno],acz[bno], + bcx[bno],bcy[bno],bcz[bno], + b,dinfo,sinfo); + }); + } else if (has_alpha && !has_beta) { + auto const& acx = m_acoefs[amrlev][mglev][0]->const_arrays(); + auto const& acy = m_acoefs[amrlev][mglev][1]->const_arrays(); + auto const& acz = m_acoefs[amrlev][mglev][2]->const_arrays(); + auto b = m_beta; + ParallelFor(nmf, [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k) + { + Array4 empty; + mlcurlcurl_gs4_alpha(i,j,k,ex[bno],ey[bno],ez[bno], + rhsx[bno],rhsy[bno],rhsz[bno], + dxinv,color, + acx[bno],acy[bno],acz[bno], + empty,empty,empty, + b,dinfo,sinfo); + }); } else { + // This branch covers scalar alpha and variable beta. + // If LU is used, we will build local solvers. + AMREX_ASSERT(!has_alpha && has_beta); auto const& bcx = m_bcoefs[amrlev][mglev][0]->const_arrays(); auto const& bcy = m_bcoefs[amrlev][mglev][1]->const_arrays(); auto const& bcz = m_bcoefs[amrlev][mglev][2]->const_arrays(); - if (m_use_pcg) { + if (use_pcg) { ParallelFor(nmf, [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k) { - mlcurlcurl_gs4(i,j,k,ex[bno],ey[bno],ez[bno], rhsx[bno],rhsy[bno],rhsz[bno], adxinv,color,bcx[bno],bcy[bno],bcz[bno], @@ -496,7 +730,6 @@ void MLCurlCurl::smooth4 (int amrlev, int mglev, MF& sol, MF const& rhs, } else { ParallelFor(nmf, [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k) { - mlcurlcurl_gs4(i,j,k,ex[bno],ey[bno],ez[bno], rhsx[bno],rhsy[bno],rhsz[bno], adxinv,color,bcx[bno],bcy[bno],bcz[bno], @@ -577,7 +810,8 @@ void MLCurlCurl::compresid (int amrlev, int mglev, MF& resid, MF const& b) const void MLCurlCurl::update_lusolver () { #if (AMREX_SPACEDIM > 1) - if (m_bcoefs[0][0][0] == nullptr) { + // There is no global LU Solver that can be built for variable alpha or beta. + if (m_bcoefs[0][0][0] == nullptr && m_acoefs[0][0][0] == nullptr) { for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) { for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev) { auto const& dxinv = this->m_geom[amrlev][mglev].InvCellSizeArray(); @@ -685,6 +919,28 @@ void MLCurlCurl::set_curvilinear_domain_bc () #endif } +void MLCurlCurl::clearAlphaMultiFab () +{ + for (auto& amrvec : m_acoefs) { + for (auto& arr : amrvec) { + for (auto& mf : arr) { + mf.reset(); + } + } + } +} + +void MLCurlCurl::clearBetaMultiFab () +{ + for (auto& amrvec : m_bcoefs) { + for (auto& arr : amrvec) { + for (auto& mf : arr) { + mf.reset(); + } + } + } +} + Real MLCurlCurl::xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const { diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H index 5658b1f8c2..c67c45b3d1 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H @@ -504,6 +504,211 @@ void mlcurlcurl_adotx_z (int i, int j, int k, Array4 const& Az, Az(i,j,k) = ccez + beta * ez(i,j,k); } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_adotx_x_alpha (int i, int j, int k, Array4 const& Ax, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Array4 const& afy, + Array4 const& afz, + Real beta, GpuArray const& dxinv) +{ +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(ey,ez,afy,afz,dxinv); + Ax(i,j,k) = beta * ex(i,j,k); +#elif (AMREX_SPACEDIM == 2) + amrex::ignore_unused(ez,afy); + + Real const alpha_hi = afz(i,j,k); + Real const alpha_lo = afz(i,j-1,k); + + Real const curlz_hi = dxinv[0] * (ey(i+1,j ,k) - ey(i ,j ,k)) + - dxinv[1] * (ex(i ,j+1,k) - ex(i ,j ,k)); + Real const curlz_lo = dxinv[0] * (ey(i+1,j-1,k) - ey(i ,j-1,k)) + - dxinv[1] * (ex(i ,j ,k) - ex(i ,j-1,k)); + + Real const ccex = dxinv[1] * (alpha_hi * curlz_hi - alpha_lo * curlz_lo); + + Ax(i,j,k) = ccex + beta * ex(i,j,k); +#else + Real const ay_lo = afy(i,j ,k-1); + Real const ay_hi = afy(i,j ,k ); + Real const az_lo = afz(i,j-1,k ); + Real const az_hi = afz(i,j ,k ); + + Real const curlz_hi = dxinv[0] * (ey(i+1,j ,k ) - ey(i ,j ,k )) + - dxinv[1] * (ex(i ,j+1,k ) - ex(i ,j ,k )); + Real const curlz_lo = dxinv[0] * (ey(i+1,j-1,k ) - ey(i ,j-1,k )) + - dxinv[1] * (ex(i ,j ,k ) - ex(i ,j-1,k )); + + Real const curly_hi = dxinv[2] * (ex(i ,j ,k+1) - ex(i ,j ,k )) + - dxinv[0] * (ez(i+1,j ,k ) - ez(i ,j ,k )); + Real const curly_lo = dxinv[2] * (ex(i ,j ,k ) - ex(i ,j ,k-1)) + - dxinv[0] * (ez(i+1,j ,k-1) - ez(i ,j ,k-1)); + + Real const ccex = dxinv[1] * (az_hi * curlz_hi - az_lo * curlz_lo) + - dxinv[2] * (ay_hi * curly_hi - ay_lo * curly_lo); + + Ax(i,j,k) = ccex + beta * ex(i,j,k); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_adotx_y_alpha (int i, int j, int k, Array4 const& Ax, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Array4 const& afx, + Array4 const& afz, + Real beta, GpuArray const& dxinv +#if (AMREX_SPACEDIM < 3) + , int coord +#endif + ) +{ +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(ex,ez,afx); + Real const dxinv0 = dxinv[0]; + Real const alpha_hi = afz(i ,j,k); + Real const alpha_lo = afz(i-1,j,k); + + // The discrezation for curvilinear coordinates does not reduce to the + // same form used by the scalar alpha case. But I think either form is fine. + + if (coord == 0) { + Real const curl_lo = dxinv0 * (ey(i ,j,k) - ey(i-1,j,k)); + Real const curl_hi = dxinv0 * (ey(i+1,j,k) - ey(i ,j,k)); + Real const ccey = dxinv0 * (alpha_lo * curl_lo - alpha_hi * curl_hi); + Ax(i,j,k) = ccey + beta * ey(i,j,k); + } else if (coord == 1) { + Real const curl_lo = dxinv0 * (Real(i)*ey(i ,j,k) - Real(i-1)*ey(i-1,j,k)) + / (Real(i)-Real(0.5)); + Real const curl_hi = dxinv0 * (Real(i+1)*ey(i+1,j,k) - Real(i)*ey(i ,j,k)) + / (Real(i)+Real(0.5)); + Real const ccey = dxinv0 * (alpha_lo * curl_lo - alpha_hi * curl_hi); + Ax(i,j,k) = ccey + beta * ey(i,j,k); + } else { + Real const curl_lo = dxinv0 * (Real(i)*ey(i ,j,k) - Real(i-1)*ey(i-1,j,k)); + Real const curl_hi = dxinv0 * (Real(i+1)*ey(i+1,j,k) - Real(i)*ey(i ,j,k)); + Real const ccey = dxinv0 * (alpha_lo*curl_lo - alpha_hi*curl_hi) / Real(i); + Ax(i,j,k) = ccey + beta * ey(i,j,k); + } +#elif (AMREX_SPACEDIM == 2) + amrex::ignore_unused(ez,afx,coord); + + Real const alpha_rt = afz(i ,j,k); + Real const alpha_lf = afz(i-1,j,k); + + Real const curlz_rt = dxinv[0] * (ey(i+1,j ,k) - ey(i ,j ,k)) + - dxinv[1] * (ex(i ,j+1,k) - ex(i ,j ,k)); + Real const curlz_lf = dxinv[0] * (ey(i ,j ,k) - ey(i-1,j ,k)) + - dxinv[1] * (ex(i-1,j+1,k) - ex(i-1,j ,k)); + + Real const ccey = dxinv[0] * (alpha_lf * curlz_lf - alpha_rt * curlz_rt); + + Ax(i,j,k) = ccey + beta * ey(i,j,k); +#else + Real const ax_lo = afx(i ,j,k-1); + Real const ax_hi = afx(i ,j,k ); + Real const az_lo = afz(i-1,j,k ); + Real const az_hi = afz(i ,j,k ); + + Real const curlz_hi = dxinv[0] * (ey(i+1,j ,k ) - ey(i ,j ,k )) + - dxinv[1] * (ex(i ,j+1,k ) - ex(i ,j ,k )); + Real const curlz_lo = dxinv[0] * (ey(i ,j ,k ) - ey(i-1,j ,k )) + - dxinv[1] * (ex(i-1,j+1,k ) - ex(i-1,j ,k )); + + Real const curlx_hi = dxinv[1] * (ez(i ,j+1,k ) - ez(i ,j ,k )) + - dxinv[2] * (ey(i ,j ,k+1) - ey(i ,j ,k )); + Real const curlx_lo = dxinv[1] * (ez(i ,j+1,k-1) - ez(i ,j ,k-1)) + - dxinv[2] * (ey(i ,j ,k ) - ey(i ,j ,k-1)); + + Real const ccey = dxinv[2] * (ax_hi * curlx_hi - ax_lo * curlx_lo) + - dxinv[0] * (az_hi * curlz_hi - az_lo * curlz_lo); + + Ax(i,j,k) = ccey + beta * ey(i,j,k); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_adotx_z_alpha (int i, int j, int k, Array4 const& Ax, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Array4 const& afx, + Array4 const& afy, + Real beta, GpuArray const& dxinv +#if (AMREX_SPACEDIM < 3) + , int coord +#endif + ) +{ +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(ex,ey,afx); + Real const dxinv0 = dxinv[0]; + Real const alpha_hi = afy(i ,j,k); + Real const alpha_lo = afy(i-1,j,k); + + // The discrezation for curvilinear coordinates does not reduce to the + // same form used by the scalar alpha case. But I think either form is fine. + + if (coord == 0) { + Real const curl_lo = dxinv0 * (ez(i ,j,k) - ez(i-1,j,k)); + Real const curl_hi = dxinv0 * (ez(i+1,j,k) - ez(i ,j,k)); + Real const ccez = dxinv0 * (alpha_lo * curl_lo - alpha_hi * curl_hi); + Ax(i,j,k) = ccez + beta * ez(i,j,k); + } else if (coord == 1) { + Real ccez; + if (i == 0) { + ccez = Real(4.0) * (ez(0,0,0) - ez(1,0,0)) * alpha_hi * dxinv0*dxinv0; + } else { + Real const curl_lo = dxinv0 * (ez(i ,j,k) - ez(i-1,j,k)); + Real const curl_hi = dxinv0 * (ez(i+1,j,k) - ez(i ,j,k)); + ccez = dxinv0 * ((Real(i)-Real(0.5))*alpha_lo*curl_lo - + (Real(i)+Real(0.5))*alpha_hi*curl_hi)/Real(i); + } + Ax(i,j,k) = ccez + beta * ez(i,j,k); + } else { + Real const curl_lo = dxinv0 * (Real(i)*ez(i ,j,k) - Real(i-1)*ez(i-1,j,k)); + Real const curl_hi = dxinv0 * (Real(i+1)*ez(i+1,j,k) - Real(i)*ez(i ,j,k)); + Real const ccez = dxinv0 * (alpha_lo*curl_lo - alpha_hi*curl_hi) / Real(i); + Ax(i,j,k) = ccez + beta * ez(i,j,k); + } +#elif (AMREX_SPACEDIM == 2) + amrex::ignore_unused(ex,ey,coord); + + Real const cx_up = dxinv[1] * (ez(i ,j+1,k) - ez(i ,j ,k)); + Real const cx_dn = dxinv[1] * (ez(i ,j ,k) - ez(i ,j-1,k)); + Real const cy_rt = -dxinv[0] * (ez(i+1,j ,k) - ez(i ,j ,k)); + Real const cy_lf = -dxinv[0] * (ez(i ,j ,k) - ez(i-1,j ,k)); + + Real const ccez = dxinv[0] * (afy(i ,j ,k) * cy_rt - afy(i-1,j ,k) * cy_lf) + - dxinv[1] * (afx(i ,j ,k) * cx_up - afx(i ,j-1,k) * cx_dn); + + Ax(i,j,k) = ccez + beta * ez(i,j,k); +#else + Real const ay_lo = afy(i-1,j ,k); + Real const ay_hi = afy(i ,j ,k); + Real const ax_lo = afx(i ,j-1,k); + Real const ax_hi = afx(i ,j ,k); + + Real const curlx_hi = dxinv[1] * (ez(i ,j+1,k ) - ez(i ,j ,k )) + - dxinv[2] * (ey(i ,j ,k+1) - ey(i ,j ,k )); + Real const curlx_lo = dxinv[1] * (ez(i ,j ,k ) - ez(i ,j-1,k )) + - dxinv[2] * (ey(i ,j-1,k+1) - ey(i ,j-1,k )); + + Real const curly_hi = dxinv[2] * (ex(i ,j ,k+1) - ex(i ,j ,k )) + - dxinv[0] * (ez(i+1,j ,k ) - ez(i ,j ,k )); + Real const curly_lo = dxinv[2] * (ex(i-1,j ,k+1) - ex(i-1,j ,k )) + - dxinv[0] * (ez(i ,j ,k ) - ez(i-1,j ,k )); + + Real const ccez = dxinv[0] * (ay_hi * curly_hi - ay_lo * curly_lo) + - dxinv[1] * (ax_hi * curlx_hi - ax_lo * curlx_lo); + + Ax(i,j,k) = ccez + beta * ez(i,j,k); +#endif +} + #if (AMREX_SPACEDIM == 1) AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_smooth_1d (int i, int j, int k, @@ -675,6 +880,226 @@ void mlcurlcurl_smooth_1d (int i, int j, int k, } } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_smooth_1d_alpha_beta (int i, int j, int k, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Array4 const& rhsx, + Array4 const& rhsy, + Array4 const& rhsz, + Array4 const& betax, + Array4 const& betay, + Array4 const& betaz, + GpuArray const& adxinv, + int color, + CurlCurlDirichletInfo const& dinfo, bool valid_x, + int coord, + Array4 const& alphay, + Array4 const& alphaz) +{ + Real dxx = adxinv[0] * adxinv[0]; + + int my_color = i%2; + if (my_color != color) { return; } + + if (valid_x) { + ex(i,j,k) = rhsx(i,j,k) / betax(i,j,k); + } + + auto cartesian_update = [&] (Array4 const& comp, + Array4 const& rhs, + Array4 const& beta, + Array4 const& alpha) + { + if (dinfo.is_dirichlet_node(i,j,k)) { return; } + Real alpha_lo = alpha(i-1,j,k); + Real alpha_hi = alpha(i ,j,k); + Real cm = -dxx * alpha_lo; + Real cp = -dxx * alpha_hi; + Real gamma = dxx * (alpha_lo + alpha_hi) + beta(i,j,k); + Real cce = cm * comp(i-1,j,k) + cp * comp(i+1,j,k); + Real res = rhs(i,j,k) - (gamma * comp(i,j,k) + cce); + comp(i,j,k) += res / gamma; + }; + + auto cylindrical_y = [&] () + { + if ((i > 0) && !dinfo.is_dirichlet_node(i,j,k)) { + Real ri = Real(i); + Real rimh = ri - Real(0.5); + Real riph = ri + Real(0.5); + Real alpha_lo = alphaz(i-1,j,k); + Real alpha_hi = alphaz(i ,j,k); + Real cm = -dxx * alpha_lo * (ri-Real(1.0)) / rimh; + Real cp = -dxx * alpha_hi * (ri+Real(1.0)) / riph; + Real gamma = dxx * (alpha_lo * ri/rimh + alpha_hi * ri/riph) + + betay(i,j,k); + Real ccey = cm * ey(i-1,j,k) + cp * ey(i+1,j,k); + Real res_y = rhsy(i,j,k) - (gamma * ey(i,j,k) + ccey); + ey(i,j,k) += res_y / gamma; + } + }; + + auto cylindrical_z = [&] () + { + if ((i > 0) && !dinfo.is_dirichlet_node(i,j,k)) { + Real ri = Real(i); + Real alpha_lo = alphay(i-1,j,k); + Real alpha_hi = alphay(i ,j,k); + Real cm = -(dxx/ri) * (ri-Real(0.5)) * alpha_lo; + Real cp = -(dxx/ri) * (ri+Real(0.5)) * alpha_hi; + Real gamma = (dxx/ri) * ((ri-Real(0.5))*alpha_lo + + (ri+Real(0.5))*alpha_hi) + + betaz(i,j,k); + Real ccez = cm * ez(i-1,j,k) + cp * ez(i+1,j,k); + Real res_z = rhsz(i,j,k) - (gamma * ez(i,j,k) + ccez); + ez(i,j,k) += res_z / gamma; + } + if ((i == 0) || (i == 1)) { + Real alpha_axis = alphay(0,j,k); + Real gamma = betaz(0,j,k) + Real(4.0)*dxx*alpha_axis; + Real num = rhsz(0,j,k) + Real(4.0)*dxx*alpha_axis * ez(1,j,k); + ez(0,j,k) = num / gamma; + } + }; + + auto spherical_update = [&] (Array4 const& comp, + Array4 const& rhs, + Array4 const& beta, + Array4 const& alpha) + { + if (i == 0 || dinfo.is_dirichlet_node(i,j,k)) { return; } + Real ri = Real(i); + Real alpha_lo = alpha(i-1,j,k); + Real alpha_hi = alpha(i ,j,k); + Real cm = -dxx * alpha_lo * Real(i-1) / ri; + Real cp = -dxx * alpha_hi * Real(i+1) / ri; + Real gamma = dxx * (alpha_lo + alpha_hi) + beta(i,j,k); + Real cce = cm * comp(i-1,j,k) + cp * comp(i+1,j,k); + Real res = rhs(i,j,k) - (gamma * comp(i,j,k) + cce); + comp(i,j,k) += res / gamma; + }; + + if (coord == 0) { + cartesian_update(ey, rhsy, betay, alphaz); + cartesian_update(ez, rhsz, betaz, alphay); + } else if (coord == 1) { + cylindrical_y(); + cylindrical_z(); + } else { + spherical_update(ey, rhsy, betay, alphaz); + spherical_update(ez, rhsz, betaz, alphay); + } +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_smooth_1d_alpha (int i, int j, int k, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Array4 const& rhsx, + Array4 const& rhsy, + Array4 const& rhsz, + Real beta, + GpuArray const& adxinv, + int color, + CurlCurlDirichletInfo const& dinfo, bool valid_x, + int coord, + Array4 const& alphay, + Array4 const& alphaz) +{ + Real dxx = adxinv[0] * adxinv[0]; + + int my_color = i%2; + if (my_color != color) { return; } + + if (valid_x) { + ex(i,j,k) = rhsx(i,j,k) / beta; + } + + auto cartesian_update = [&] (Array4 const& comp, + Array4 const& rhs, + Array4 const& alpha) + { + if (dinfo.is_dirichlet_node(i,j,k)) { return; } + Real alpha_lo = alpha(i-1,j,k); + Real alpha_hi = alpha(i ,j,k); + Real cm = -dxx * alpha_lo; + Real cp = -dxx * alpha_hi; + Real gamma = dxx * (alpha_lo + alpha_hi) + beta; + Real cce = cm * comp(i-1,j,k) + cp * comp(i+1,j,k); + Real res = rhs(i,j,k) - (gamma * comp(i,j,k) + cce); + comp(i,j,k) += res / gamma; + }; + + auto cylindrical_y = [&] () + { + if ((i > 0) && !dinfo.is_dirichlet_node(i,j,k)) { + Real ri = Real(i); + Real rimh = ri - Real(0.5); + Real riph = ri + Real(0.5); + Real alpha_lo = alphaz(i-1,j,k); + Real alpha_hi = alphaz(i ,j,k); + Real cm = -dxx * alpha_lo * (ri-Real(1.0)) / rimh; + Real cp = -dxx * alpha_hi * (ri+Real(1.0)) / riph; + Real gamma = dxx * (alpha_lo * ri/rimh + alpha_hi * ri/riph) + beta; + Real ccey = cm * ey(i-1,j,k) + cp * ey(i+1,j,k); + Real res_y = rhsy(i,j,k) - (gamma * ey(i,j,k) + ccey); + ey(i,j,k) += res_y / gamma; + } + }; + + auto cylindrical_z = [&] () + { + if ((i > 0) && !dinfo.is_dirichlet_node(i,j,k)) { + Real ri = Real(i); + Real alpha_lo = alphay(i-1,j,k); + Real alpha_hi = alphay(i ,j,k); + Real cm = -(dxx/ri) * (ri-Real(0.5)) * alpha_lo; + Real cp = -(dxx/ri) * (ri+Real(0.5)) * alpha_hi; + Real gamma = (dxx/ri) * ((ri-Real(0.5))*alpha_lo + + (ri+Real(0.5))*alpha_hi) + beta; + Real ccez = cm * ez(i-1,j,k) + cp * ez(i+1,j,k); + Real res_z = rhsz(i,j,k) - (gamma * ez(i,j,k) + ccez); + ez(i,j,k) += res_z / gamma; + } + if ((i == 0) || (i == 1)) { + Real alpha_axis = alphay(0,j,k); + Real gamma = beta + Real(4.0)*dxx*alpha_axis; + Real num = rhsz(0,j,k) + Real(4.0)*dxx*alpha_axis * ez(1,j,k); + ez(0,j,k) = num / gamma; + } + }; + + auto spherical_update = [&] (Array4 const& comp, + Array4 const& rhs, + Array4 const& alpha) + { + if (i == 0 || dinfo.is_dirichlet_node(i,j,k)) { return; } + Real ri = Real(i); + Real alpha_lo = alpha(i-1,j,k); + Real alpha_hi = alpha(i ,j,k); + Real cm = -dxx * alpha_lo * Real(i-1) / ri; + Real cp = -dxx * alpha_hi * Real(i+1) / ri; + Real gamma = dxx * (alpha_lo + alpha_hi) + beta; + Real cce = cm * comp(i-1,j,k) + cp * comp(i+1,j,k); + Real res = rhs(i,j,k) - (gamma * comp(i,j,k) + cce); + comp(i,j,k) += res / gamma; + }; + + if (coord == 0) { + cartesian_update(ey, rhsy, alphaz); + cartesian_update(ez, rhsz, alphay); + } else if (coord == 1) { + cylindrical_y(); + cylindrical_z(); + } else { + spherical_update(ey, rhsy, alphaz); + spherical_update(ez, rhsz, alphay); + } +} + #endif @@ -725,7 +1150,7 @@ void mlcurlcurl_gs4_lu (int i, int j, int k, Real dxy = adxinv[0]*adxinv[1]; - GpuArray b + GpuArray b {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) + ex(i-1,j+1,k )) + dxy * ( ey(i-1,j-1,k ) @@ -898,7 +1323,7 @@ void mlcurlcurl_gs4 (int i, int j, int k, Real dxy = adxinv[0]*adxinv[1]; - GpuArray b + GpuArray b {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) + ex(i-1,j+1,k )) + dxy * ( ey(i-1,j-1,k ) @@ -1178,6 +1603,321 @@ void mlcurlcurl_gs4 (int i, int j, int k, #endif } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void mlcurlcurl_gs4_alpha (int i, int j, int k, + Array4 const& ex, + Array4 const& ey, + Array4 const& ez, + Array4 const& rhsx, + Array4 const& rhsy, + Array4 const& rhsz, + GpuArray const& dxinv, + int color, + Array4 const& alphax, + Array4 const& alphay, + Array4 const& alphaz, + Array4 const& betax, + Array4 const& betay, + Array4 const& betaz, Real bscalar, + CurlCurlDirichletInfo const& dinfo, + CurlCurlSymmetryInfo const& sinfo) +{ + if (dinfo.is_dirichlet_node(i,j,k)) { return; } + + int my_color = i%2 + 2*(j%2); + if (k%2 != 0) { + my_color = 3 - my_color; + } + +#if (AMREX_SPACEDIM == 2) + + Real dxx = dxinv[0] * dxinv[0]; + Real dyy = dxinv[1] * dxinv[1]; + + if (((my_color == 0 || my_color == 3) && (color == 0 || color == 3)) || + ((my_color == 1 || my_color == 2) && (color == 1 || color == 2))) + { + Real beta = betaz ? betaz(i,j,k) : bscalar; + Real gamma = dxx*(alphay(i-1,j ,k)+alphay(i,j,k)) + + dyy*(alphax(i ,j-1,k)+alphax(i,j,k)) + beta; + Real ccez = - dxx * (ez(i-1,j ,k ) * alphay(i-1,j ,k) + + ez(i+1,j ,k ) * alphay(i ,j ,k)) + - dyy * (ez(i ,j-1,k ) * alphax(i ,j-1,k) + + ez(i ,j+1,k ) * alphax(i ,j ,k)); + Real res = rhsz(i,j,k) - (gamma*ez(i,j,k) + ccez); + constexpr Real omega = Real(1.15); + ez(i,j,k) += omega/gamma * res; + } + + if (my_color != color) { return; } + + Real dxy = dxinv[0]*dxinv[1]; + + Real amm = alphaz(i-1,j-1,k); + Real am0 = alphaz(i-1,j ,k); + Real a0m = alphaz(i ,j-1,k); + Real a00 = alphaz(i ,j ,k); + + GpuArray b + {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) * amm + + ex(i-1,j+1,k ) * am0 ) + + dxy * ( ey(i-1,j-1,k ) * amm + -ey(i-1,j ,k ) * am0 )), + rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) * a0m + + ex(i ,j+1,k ) * a00 ) + + dxy * (-ey(i+1,j-1,k ) * a0m + +ey(i+1,j ,k ) * a00 )), + rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) * amm + + ey(i+1,j-1,k ) * a0m ) + + dxy * ( ex(i-1,j-1,k ) * amm + -ex(i ,j-1,k ) * a0m )), + rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) * am0 + + ey(i+1,j ,k ) * a00 ) + + dxy * (-ex(i-1,j+1,k ) * am0 + +ex(i ,j+1,k ) * a00))}; + + GpuArray beta; + + Real bxm = bscalar; + Real bx0 = bscalar; + if (betax) { + bxm = betax(i-1,j,k); + bx0 = betax(i ,j,k); + } + Real bym = bscalar; + Real by0 = bscalar; + if (betay) { + bym = betay(i,j-1,k); + by0 = betay(i,j ,k); + } + + if (sinfo.xlo_is_symmetric(i)) { + b[0] = -b[1]; + beta[0] = beta[1] = bx0; + } else if (sinfo.xhi_is_symmetric(i)) { + b[1] = -b[0]; + beta[0] = beta[1] = bxm; + } else { + beta[0] = bxm; + beta[1] = bx0; + } + + if (sinfo.ylo_is_symmetric(j)) { + b[2] = -b[3]; + beta[2] = beta[3] = by0; + } else if (sinfo.yhi_is_symmetric(j)) { + b[3] = -b[2]; + beta[2] = beta[3] = bym; + } else { + beta[2] = bym; + beta[3] = by0; + } + + Real diagInv[4] = {Real(1.0) / (dyy*(amm+am0) + beta[0]), + Real(1.0) / (dyy*(a0m+a00) + beta[1]), + Real(1.0) / (dxx*(amm+a0m) + beta[2]), + Real(1.0) / (dxx*(am0+a00) + beta[3])}; + auto precond = [&] (Real * AMREX_RESTRICT z, + Real const* AMREX_RESTRICT r) + { + for (int m = 0; m < 4; ++m) { z[m] = r[m] * diagInv[m]; } + }; + auto mat = [&] (Real * AMREX_RESTRICT Av, + Real const* AMREX_RESTRICT v) + { + Av[0] = (dyy*(am0+amm) + beta[0]) * v[0] - dxy*amm * v[2] + dxy*am0 * v[3]; + Av[1] = (dyy*(a00+a0m) + beta[1]) * v[1] + dxy*a0m * v[2] - dxy*a00 * v[3]; + Av[2] = -dxy*amm * v[0] + dxy*a0m * v[1] + (dxx*(amm+a0m) + beta[2]) * v[2]; + Av[3] = dxy*am0 * v[0] - dxy*a00 * v[1] + (dxx*(am0+a00) + beta[3]) * v[3]; + }; + Real sol[4] = {0, 0, 0, 0}; + pcg_solve<4>(sol, b.data(), mat, precond, 8, Real(1.e-8)); + ex(i-1,j ,k ) = sol[0]; + ex(i ,j ,k ) = sol[1]; + ey(i ,j-1,k ) = sol[2]; + ey(i ,j ,k ) = sol[3]; + +#elif (AMREX_SPACEDIM == 3) + + if (my_color != color) { return; } + + Real dxx = dxinv[0]*dxinv[0]; + 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 dyz = dxinv[1]*dxinv[2]; + + Real const az_mm = alphaz(i-1,j-1,k); + Real const az_m0 = alphaz(i-1,j ,k); + Real const az_0m = alphaz(i ,j-1,k); + Real const az_00 = alphaz(i ,j ,k); + + Real const ay_mm = alphay(i-1,j ,k-1); + Real const ay_m0 = alphay(i-1,j ,k ); + Real const ay_0m = alphay(i ,j ,k-1); + Real const ay_00 = alphay(i ,j ,k ); + + Real const ax_mm = alphax(i ,j-1,k-1); + Real const ax_m0 = alphax(i ,j-1,k ); + Real const ax_0m = alphax(i ,j ,k-1); + Real const ax_00 = alphax(i ,j ,k ); + + GpuArray b + {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) * az_mm + + ex(i-1,j+1,k ) * az_m0) + - dzz * ( ex(i-1,j ,k-1) * ay_mm + + ex(i-1,j ,k+1) * ay_m0) + + dxy * ( ey(i-1,j-1,k ) * az_mm + -ey(i-1,j ,k ) * az_m0) + + dxz * ( ez(i-1,j ,k-1) * ay_mm + -ez(i-1,j ,k ) * ay_m0)), + rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) * az_0m + + ex(i ,j+1,k ) * az_00) + - dzz * ( ex(i ,j ,k-1) * ay_0m + + ex(i ,j ,k+1) * ay_00) + + dxy * (-ey(i+1,j-1,k ) * az_0m + +ey(i+1,j ,k ) * az_00) + + dxz * (-ez(i+1,j ,k-1) * ay_0m + +ez(i+1,j ,k ) * ay_00)), + rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) * az_mm + + ey(i+1,j-1,k ) * az_0m) + - dzz * ( ey(i ,j-1,k-1) * ax_mm + + ey(i ,j-1,k+1) * ax_m0) + + dxy * ( ex(i-1,j-1,k ) * az_mm + -ex(i ,j-1,k ) * az_0m) + + dyz * ( ez(i ,j-1,k-1) * ax_mm + -ez(i ,j-1,k ) * ax_m0)), + rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) * az_m0 + + ey(i+1,j ,k ) * az_00) + - dzz * ( ey(i ,j ,k-1) * ax_0m + + ey(i ,j ,k+1) * ax_00) + + dxy * (-ex(i-1,j+1,k ) * az_m0 + +ex(i ,j+1,k ) * az_00) + + dyz * (-ez(i ,j+1,k-1) * ax_0m + +ez(i ,j+1,k ) * ax_00)), + rhsz(i,j,k-1) - (-dxx * ( ez(i-1,j ,k-1) * ay_mm + + ez(i+1,j ,k-1) * ay_0m) + - dyy * ( ez(i ,j-1,k-1) * ax_mm + + ez(i ,j+1,k-1) * ax_0m) + + dxz * ( ex(i-1,j ,k-1) * ay_mm + -ex(i ,j ,k-1) * ay_0m) + + dyz * ( ey(i ,j-1,k-1) * ax_mm + -ey(i ,j ,k-1) * ax_0m)), + rhsz(i,j,k ) - (-dxx * ( ez(i-1,j ,k ) * ay_m0 + + ez(i+1,j ,k ) * ay_00) + - dyy * ( ez(i ,j-1,k ) * ax_m0 + + ez(i ,j+1,k ) * ax_00) + + dxz * (-ex(i-1,j ,k+1) * ay_m0 + +ex(i ,j ,k+1) * ay_00) + + dyz * (-ey(i ,j-1,k+1) * ax_m0 + +ey(i ,j ,k+1) * ax_00))}; + + auto get_bx = [&] (int ii, int jj, int kk) -> Real + { + return betax ? betax(ii,jj,kk) : bscalar; + }; + auto get_by = [&] (int ii, int jj, int kk) -> Real + { + return betay ? betay(ii,jj,kk) : bscalar; + }; + auto get_bz = [&] (int ii, int jj, int kk) -> Real + { + return betaz ? betaz(ii,jj,kk) : bscalar; + }; + + GpuArray beta; + + if (sinfo.xlo_is_symmetric(i)) { + b[0] = -b[1]; + Real bx = get_bx(i,j,k); + beta[0] = beta[1] = bx; + } else if (sinfo.xhi_is_symmetric(i)) { + b[1] = -b[0]; + Real bx = get_bx(i-1,j,k); + beta[0] = beta[1] = bx; + } else { + beta[0] = get_bx(i-1,j,k); + beta[1] = get_bx(i ,j,k); + } + + if (sinfo.ylo_is_symmetric(j)) { + b[2] = -b[3]; + Real by = get_by(i,j,k); + beta[2] = beta[3] = by; + } else if (sinfo.yhi_is_symmetric(j)) { + b[3] = -b[2]; + Real by = get_by(i,j-1,k); + beta[2] = beta[3] = by; + } else { + beta[2] = get_by(i,j-1,k); + beta[3] = get_by(i,j ,k); + } + + if (sinfo.zlo_is_symmetric(k)) { + b[4] = -b[5]; + Real bz = get_bz(i,j,k); + beta[4] = beta[5] = bz; + } else if (sinfo.zhi_is_symmetric(k)) { + b[5] = -b[4]; + Real bz = get_bz(i,j,k-1); + beta[4] = beta[5] = bz; + } else { + beta[4] = get_bz(i,j,k-1); + beta[5] = get_bz(i,j,k ); + } + + Real diag0 = dyy*(az_mm+az_m0) + dzz*(ay_mm+ay_m0) + beta[0]; + Real diag1 = dyy*(az_0m+az_00) + dzz*(ay_0m+ay_00) + beta[1]; + Real diag2 = dxx*(az_mm+az_0m) + dzz*(ax_mm+ax_m0) + beta[2]; + Real diag3 = dxx*(az_m0+az_00) + dzz*(ax_0m+ax_00) + beta[3]; + Real diag4 = dxx*(ay_mm+ay_0m) + dyy*(ax_mm+ax_0m) + beta[4]; + Real diag5 = dxx*(ay_m0+ay_00) + dyy*(ax_m0+ax_00) + beta[5]; + + Real diagInv[6] = {Real(1.0) / diag0, + Real(1.0) / diag1, + Real(1.0) / diag2, + Real(1.0) / diag3, + Real(1.0) / diag4, + Real(1.0) / diag5}; + auto precond = [&] (Real * AMREX_RESTRICT z, + Real const* AMREX_RESTRICT r) + { + for (int m = 0; m < 6; ++m) { z[m] = r[m] * diagInv[m]; } + }; + auto mat = [&] (Real * AMREX_RESTRICT Av, + Real const* AMREX_RESTRICT v) + { + Av[0] = diag0 * v[0] + - dxy * az_mm * v[2] + dxy * az_m0 * v[3] + - dxz * ay_mm * v[4] + dxz * ay_m0 * v[5]; + Av[1] = diag1 * v[1] + + dxy * az_0m * v[2] - dxy * az_00 * v[3] + + dxz * ay_0m * v[4] - dxz * ay_00 * v[5]; + Av[2] = diag2 * v[2] + - dxy * az_mm * v[0] + dxy * az_0m * v[1] + - dyz * ax_mm * v[4] + dyz * ax_m0 * v[5]; + Av[3] = diag3 * v[3] + + dxy * az_m0 * v[0] - dxy * az_00 * v[1] + + dyz * ax_0m * v[4] - dyz * ax_00 * v[5]; + Av[4] = diag4 * v[4] + - dxz * ay_mm * v[0] + dxz * ay_0m * v[1] + - dyz * ax_mm * v[2] + dyz * ax_0m * v[3]; + Av[5] = diag5 * v[5] + + dxz * ay_m0 * v[0] - dxz * ay_00 * v[1] + + dyz * ax_m0 * v[2] - dyz * ax_00 * v[3]; + }; + Real sol[6] = {0, 0, 0, 0, 0, 0}; + pcg_solve<6>(sol, b.data(), mat, precond, 8, Real(1.e-8)); + ex(i-1,j ,k ) = sol[0]; + ex(i ,j ,k ) = sol[1]; + ey(i ,j-1,k ) = sol[2]; + ey(i ,j ,k ) = sol[3]; + ez(i ,j ,k-1) = sol[4]; + ez(i ,j ,k ) = sol[5]; +#endif +} + #endif AMREX_GPU_DEVICE AMREX_FORCE_INLINE diff --git a/Src/LinearSolvers/MLMG/AMReX_PCGSolver.H b/Src/LinearSolvers/MLMG/AMReX_PCGSolver.H index 1b9aa55426..970a1b66cf 100644 --- a/Src/LinearSolvers/MLMG/AMReX_PCGSolver.H +++ b/Src/LinearSolvers/MLMG/AMReX_PCGSolver.H @@ -34,7 +34,7 @@ int pcg_solve (T* AMREX_RESTRICT x, T* AMREX_RESTRICT r, int iter = 0; T rho_prev = T(1.0); // initialized to quiet gcc warning - T p[N]; + T p[N] = {}; // initialized to quiet gcc warning for (iter = 1; iter <= maxiter; ++iter) { T z[N]; precond(z, r); diff --git a/Tests/LinearSolvers/CurlCurl/MyTest.H b/Tests/LinearSolvers/CurlCurl/MyTest.H index 1326894323..16aa986108 100644 --- a/Tests/LinearSolvers/CurlCurl/MyTest.H +++ b/Tests/LinearSolvers/CurlCurl/MyTest.H @@ -43,11 +43,13 @@ private: amrex::Array solution; amrex::Array exact; amrex::Array rhs; + amrex::MultiFab alpha_node; amrex::Real beta_factor = 1.e-2; amrex::Real alpha = 1.0; amrex::Real beta = 1.0; bool variable_beta = false; + bool variable_alpha = false; }; #endif diff --git a/Tests/LinearSolvers/CurlCurl/MyTest.cpp b/Tests/LinearSolvers/CurlCurl/MyTest.cpp index 2a3933c9d4..96acc98dc7 100644 --- a/Tests/LinearSolvers/CurlCurl/MyTest.cpp +++ b/Tests/LinearSolvers/CurlCurl/MyTest.cpp @@ -18,6 +18,9 @@ MyTest::MyTest () void MyTest::solve () { + Array names{"Ex", "Ey", "Ez"}; + auto dvol = AMREX_D_TERM(geom.CellSize(0), *geom.CellSize(1), *geom.CellSize(2)); + LPInfo info; info.setAgglomeration(agglomeration); info.setConsolidation(consolidation); @@ -44,6 +47,9 @@ MyTest::solve () bcoef.data()+1, bcoef.data()+2}}); } + if (variable_alpha) { + mlcc.setAlpha({&alpha_node}); + } mlcc.prepareRHS({&rhs}); if (use_pcg) { mlcc.setUsePCG(true); } @@ -53,6 +59,30 @@ MyTest::solve () mlmg.setMaxIter(max_iter); mlmg.setVerbose(verbose); mlmg.setBottomVerbose(bottom_verbose); + + { + Array ax_exact; + Array residual; + for (int idim = 0; idim < 3; ++idim) { + ax_exact[idim].define(rhs[idim].boxArray(), rhs[idim].DistributionMap(), 1, 0); + residual[idim].define(rhs[idim].boxArray(), rhs[idim].DistributionMap(), 1, 0); + MultiFab::Copy(residual[idim], rhs[idim], 0, 0, 1, 0); + } + mlmg.apply({&ax_exact}, {&exact}); + amrex::Print() << " Residual norms (b - A x_exact), expect O(dx^2):\n"; + for (int idim = 0; idim < 3; ++idim) { + MultiFab::Subtract(residual[idim], ax_exact[idim], 0, 0, 1, 0); + } + mlcc.setDirichletNodesToZero(0,0,residual); + for (int idim = 0; idim < 3; ++idim) { + auto r0 = residual[idim].norminf(); + auto r1 = residual[idim].norm1(0, geom.periodicity()) * dvol; + auto r2 = residual[idim].norm2(0, geom.periodicity()) * std::sqrt(dvol); + amrex::Print() << " " << names[idim] << ": " + << r0 << " " << r1 << " " << r2 << '\n'; + } + } + for (auto& mf : solution) { mf.setVal(Real(0)); } @@ -74,8 +104,6 @@ MyTest::solve () } amrex::Print() << " Number of cells: " << n_cell << '\n'; - 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(); @@ -117,6 +145,7 @@ MyTest::readParameters () pp.query("beta_factor", beta_factor); pp.query("alpha", alpha); pp.query("variable_beta", variable_beta); + pp.query("variable_alpha", variable_alpha); } void @@ -148,6 +177,9 @@ MyTest::initData () exact [idim].define(ba,dmap,1,1); rhs [idim].define(ba,dmap,1,0); } + if (variable_alpha) { + alpha_node.define(amrex::convert(grids,IntVect(1)), dmap,1,0); + } initProb(); diff --git a/Tests/LinearSolvers/CurlCurl/initProb.cpp b/Tests/LinearSolvers/CurlCurl/initProb.cpp index 3f7900a98c..82fff3f3dd 100644 --- a/Tests/LinearSolvers/CurlCurl/initProb.cpp +++ b/Tests/LinearSolvers/CurlCurl/initProb.cpp @@ -4,6 +4,8 @@ using namespace amrex; +enum class CoordID { Cartesian, Cyl1D, Cyl2D, Sph1D }; + void MyTest::initProb () { @@ -13,7 +15,6 @@ MyTest::initProb () const auto b = beta; const auto ndhi = geom.Domain().bigEnd()+1; - enum class CoordID { Cartesian, Cyl1D, Cyl2D, Sph1D }; auto cid = CoordID::Cartesian; #if (AMREX_SPACEDIM == 1) if (this->coord == 1) { @@ -43,18 +44,22 @@ MyTest::initProb () GpuArray,3> solfab{solution[0].array(mfi), solution[1].array(mfi), solution[2].array(mfi)}; + Array4 alphafab; + if (variable_alpha) { + alphafab = alpha_node.array(mfi); + } amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { #if (AMREX_SPACEDIM == 1) if (cid == CoordID::Sph1D) { - actual_init_prob_sph1d(i,j,k,rhsfab,solfab,prob_lo,dx,a,b,ndhi); + actual_init_prob_sph1d(i,j,k,rhsfab,solfab,prob_lo,dx,a,b,ndhi,alphafab); } else if (cid == CoordID::Cyl1D) { - actual_init_prob_cyl1d(i,j,k,rhsfab,solfab,prob_lo,dx,a,b,ndhi); + actual_init_prob_cyl1d(i,j,k,rhsfab,solfab,prob_lo,dx,a,b,ndhi,alphafab); } else #endif { - actual_init_prob(i,j,k,rhsfab,solfab,prob_lo,dx,a,b); + actual_init_prob(i,j,k,rhsfab,solfab,prob_lo,dx,a,b,alphafab); } }); } diff --git a/Tests/LinearSolvers/CurlCurl/initProb_K.H b/Tests/LinearSolvers/CurlCurl/initProb_K.H index 6399be3fee..f0a7424dfc 100644 --- a/Tests/LinearSolvers/CurlCurl/initProb_K.H +++ b/Tests/LinearSolvers/CurlCurl/initProb_K.H @@ -9,7 +9,8 @@ void actual_init_prob (int i, int j, int k, amrex::GpuArray,3> const& sol, amrex::GpuArray const& problo, amrex::GpuArray const& dx, - amrex::Real alpha, amrex::Real beta) + amrex::Real alpha, amrex::Real beta, + amrex::Array4 const& alpha_node) { using namespace amrex; @@ -26,6 +27,18 @@ void actual_init_prob (int i, int j, int k, Real zcc = znd + Real(0.5)*dx[2]; #endif + auto get_alpha = [&] (AMREX_D_DECL(Real x, Real y, Real z)) + { + Real a = AMREX_D_TERM(std::cos(Real(2.0)*pi*x), + * std::cos(Real(3.0)*pi*y), + * std::cos(Real(4.0)*pi*z)); + return Real(1.0) + Real(0.3)*a; + }; + + if (alpha_node && alpha_node.contains(i,j,k)) { + alpha_node(i,j,k) = get_alpha(AMREX_D_DECL(xnd,ynd,znd)); + } + if (sol[0].contains(i,j,k)) { Real x = xcc; Real Ex = std::sin(pi*x); @@ -69,8 +82,8 @@ void actual_init_prob (int i, int j, int k, } if (rhs[0].contains(i,j,k)) { -#if (AMREX_SPACEDIM > 1) Real x = xcc; +#if (AMREX_SPACEDIM > 1) Real y = ynd; #endif #if (AMREX_SPACEDIM == 1) @@ -85,7 +98,27 @@ void actual_init_prob (int i, int j, int k, - 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); + + Real a = (alpha_node) ? get_alpha(AMREX_D_DECL(x,y,z)) : alpha; + rhs[0](i,j,k) = a*cce + beta*sol[0](i,j,k); + +#if (AMREX_SPACEDIM >= 2) + if (alpha_node) { +#if (AMREX_SPACEDIM == 2) + Real rhs2 = Real(-3.)*pi*std::cos(Real(2.)*pi*x)*std::sin(Real(3.)*pi*y) * + (Real(-2.5)*pi*std::sin(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y) + + Real(-2.5)*pi*std::sin(pi*x)*std::cos(Real(2.5)*pi*y)); +#elif (AMREX_SPACEDIM == 3) + Real rhs2 = Real(-3.)*pi*std::cos(Real(2)*pi*x)*std::sin(Real(3)*pi*y)*std::cos(Real(4)*pi*z) * + (Real(-2.5)*pi*std::sin(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*std::sin(pi*x)*std::cos(Real(2.5)*pi*y)*std::sin(Real(2)*pi*z+Real(1./3.)*pi)) + + Real(4.)*pi*std::cos(Real(2)*pi*x)*std::cos(Real(3)*pi*y)*std::sin(Real(4)*pi*z) * + (Real(2)*pi*std::sin(pi*x)*std::sin(Real(2.5)*pi*y)*std::cos(Real(2)*pi*z+Real(1./3.)*pi) + + Real(3.5)*pi*std::sin(Real(3.5)*pi*x)*std::sin(Real(3.5)*pi*y)*std::sin(Real(4)*pi*z+Real(1./6.)*pi)); +#endif + rhs[0](i,j,k) += Real(0.3)*rhs2; + } +#endif } if (rhs[1].contains(i,j,k)) { @@ -105,7 +138,27 @@ void actual_init_prob (int i, int j, int k, + 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); + + Real a = (alpha_node) ? get_alpha(AMREX_D_DECL(x,y,z)) : alpha; + rhs[1](i,j,k) = a*cce + beta*sol[1](i,j,k); + + if (alpha_node) { +#if (AMREX_SPACEDIM == 1) + Real rhs2 = Real(-5.0)*pi*pi* std::sin(Real(2.0)*pi*x)* std::sin(Real(2.5)*pi*x); +#elif (AMREX_SPACEDIM == 2) + Real rhs2 = Real(2)*pi*std::sin(Real(2.)*pi*x)*std::cos(Real(3.)*pi*y) * + (Real(-2.5)*pi*std::sin(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y) + + Real(-2.5)*pi*std::sin(pi*x)*std::cos(Real(2.5)*pi*y)); +#elif (AMREX_SPACEDIM == 3) + Real rhs2 = Real(-4.)*pi*std::cos(Real(2)*pi*x)*std::cos(Real(3.)*pi*y)*std::sin(Real(4.)*pi*z) * + (Real(3.5)*pi*std::cos(Real(3.5)*pi*x)*std::cos(Real(3.5)*pi*y)*std::sin(Real(4.)*pi*z+Real(1./6.)*pi) + + Real(-4.)*pi*std::cos(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y)*std::cos(Real(4.)*pi*z+Real(0.25)*pi)) + + Real(2)*pi*std::sin(Real(2)*pi*x)*std::cos(Real(3.)*pi*y)*std::cos(Real(4.)*pi*z) * + (Real(-2.5)*pi*std::sin(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*std::sin(pi*x)*std::cos(Real(2.5)*pi*y)*std::sin(Real(2)*pi*z+Real(1./3.)*pi)); +#endif + rhs[1](i,j,k) += Real(0.3)*rhs2; + } } if (rhs[2].contains(i,j,k)) { @@ -123,12 +176,34 @@ void actual_init_prob (int i, int j, int k, + 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); + + Real a = (alpha_node) ? get_alpha(AMREX_D_DECL(x,y,z)) : alpha; + rhs[2](i,j,k) = a*cce + beta*sol[2](i,j,k); + + if (alpha_node) { +#if (AMREX_SPACEDIM == 1) + Real rhs2 = Real(-7.0)*pi*pi*std::sin(Real(2.0)*pi*x)*std::sin(Real(3.5)*pi*x); +#elif (AMREX_SPACEDIM == 2) + Real rhs2 = Real(-2)*pi*Real(3.5)*pi * + std::sin(Real(2.)*pi*x)*std::cos(Real(3.)*pi*y)* + std::sin(Real(3.5)*pi*x)*std::sin(Real(3.5)*pi*y) + + Real(3)*pi*Real(3.5)*pi * + std::cos(Real(2.)*pi*x)*std::sin(Real(3.)*pi*y)* + std::cos(Real(3.5)*pi*x)*std::cos(Real(3.5)*pi*y); +#elif (AMREX_SPACEDIM == 3) + Real rhs2 = Real(-2.)*pi*std::sin(Real(2.)*pi*x)*std::cos(Real(3.)*pi*y)*std::cos(Real(4.)*pi*z) * + (Real(2.)*pi*std::sin(pi*x)*std::sin(Real(2.5)*pi*y)*std::cos(Real(2.)*pi*z+Real(1./3.)*pi) + + Real(3.5)*pi*std::sin(Real(3.5)*pi*x)*std::sin(Real(3.5)*pi*y)*std::sin(Real(4.)*pi*z+Real(1./6.)*pi)) + + Real(3.)*pi*std::cos(Real(2.)*pi*x)*std::sin(Real(3.)*pi*y)*std::cos(Real(4.)*pi*z) * + (Real(3.5)*pi*std::cos(Real(3.5)*pi*x)*std::cos(Real(3.5)*pi*y)*std::sin(Real(4.)*pi*z+Real(1./6.)*pi) - + Real(4.)*pi*std::cos(Real(2.5)*pi*x)*std::sin(Real(3.)*pi*y)*std::cos(Real(4.)*pi*z+Real(0.25)*pi)); +#endif + rhs[2](i,j,k) += Real(0.3)*rhs2; + } } } #if (AMREX_SPACEDIM == 1) - AMREX_GPU_DEVICE AMREX_FORCE_INLINE void actual_init_prob_sph1d (int i, int j, int k, amrex::GpuArray,3> const& rhs, @@ -136,15 +211,30 @@ void actual_init_prob_sph1d (int i, int j, int k, amrex::GpuArray const& problo, amrex::GpuArray const& dx, amrex::Real alpha, amrex::Real beta, - amrex::IntVect const& ndhi) + amrex::IntVect const& ndhi, + amrex::Array4 const& alpha_node) { using namespace amrex; constexpr Real pi = amrex::Math::pi(); + auto alpha_val = [=] AMREX_GPU_DEVICE (Real r) noexcept + { + return Real(1.0) + Real(0.3)*std::cos(Real(2.0)*pi*r); + }; + + auto dalpha_val = [=] AMREX_GPU_DEVICE (Real r) noexcept + { + return Real(-0.6)*pi*std::sin(Real(2.0)*pi*r); + }; + Real xnd = problo[0] + Real(i)*dx[0]; Real xcc = xnd + Real(0.5)*dx[0]; + if (alpha_node && alpha_node.contains(i,j,k)) { + alpha_node(i,j,k) = alpha_val(xnd); + } + if (sol[0].contains(i,j,k)) { Real x = xcc; Real Ex = std::sin(pi*x); @@ -174,7 +264,18 @@ void actual_init_prob_sph1d (int i, int j, int k, } else { auto s = Real(2.0)*pi; Real cce = -Real(2.0)*s*std::cos(s*x)/x + s*s*std::sin(s*x); - rhs[1](i,j,k) = alpha*cce + beta*sol[1](i,j,k); + Real alpha_loc = alpha; + Real dalpha = Real(0.0); + if (alpha_node) { + alpha_loc = alpha_val(x); + dalpha = dalpha_val(x); + } + Real rhs_val = alpha_loc*cce + beta*sol[1](i,j,k); + if (alpha_node) { + // variable-alpha correction = +alpha'(r)*(curl E)_phi = -alpha'(r)*(dEy/dr + Ey/r) + rhs_val -= dalpha * (s*std::cos(s*x) + std::sin(s*x)/x); + } + rhs[1](i,j,k) = rhs_val; } } @@ -185,7 +286,18 @@ void actual_init_prob_sph1d (int i, int j, int k, } else { auto s = Real(3.0)*pi; Real cce = -Real(2.0)*s*std::cos(s*x)/x + s*s*std::sin(s*x); - rhs[2](i,j,k) = alpha*cce + beta*sol[2](i,j,k); + Real alpha_loc = alpha; + Real dalpha = Real(0.0); + if (alpha_node) { + alpha_loc = alpha_val(x); + dalpha = dalpha_val(x); + } + Real rhs_val = alpha_loc*cce + beta*sol[2](i,j,k); + if (alpha_node) { + // variable-alpha correction = +alpha'(r)*(curl E)_theta = -alpha'(r)*(dEz/dr + Ez/r) + rhs_val -= dalpha * (s*std::cos(s*x) + std::sin(s*x)/x); + } + rhs[2](i,j,k) = rhs_val; } } } @@ -197,15 +309,30 @@ void actual_init_prob_cyl1d (int i, int j, int k, amrex::GpuArray const& problo, amrex::GpuArray const& dx, amrex::Real alpha, amrex::Real beta, - amrex::IntVect const& ndhi) + amrex::IntVect const& ndhi, + amrex::Array4 const& alpha_node) { using namespace amrex; constexpr Real pi = amrex::Math::pi(); + auto alpha_val = [=] AMREX_GPU_DEVICE (Real r) noexcept + { + return Real(1.0) + Real(0.3)*std::cos(Real(2.0)*pi*r); + }; + + auto dalpha_val = [=] AMREX_GPU_DEVICE (Real r) noexcept + { + return Real(-0.6)*pi*std::sin(Real(2.0)*pi*r); + }; + Real xnd = problo[0] + Real(i)*dx[0]; Real xcc = xnd + Real(0.5)*dx[0]; + if (alpha_node && alpha_node.contains(i,j,k)) { + alpha_node(i,j,k) = alpha_val(xnd); + } + if (sol[0].contains(i,j,k)) { Real x = xcc; Real Ex = std::sin(pi*x); @@ -235,7 +362,18 @@ void actual_init_prob_cyl1d (int i, int j, int k, } else { auto s = Real(2.0)*pi; Real cce = std::sin(s*x)*(s*s+Real(1.0)/(x*x)) - std::cos(s*x)*s/x; - rhs[1](i,j,k) = alpha*cce + beta*sol[1](i,j,k); + Real alpha_loc = alpha; + Real dalpha = Real(0.0); + if (alpha_node) { + alpha_loc = alpha_val(x); + dalpha = dalpha_val(x); + } + Real rhs_val = alpha_loc*cce + beta*sol[1](i,j,k); + if (alpha_node) { + // variable-alpha correction = +alpha'(r)*(curl E)_z = -alpha'(r)*(dEy/dr + Ey/r) + rhs_val -= dalpha * (s*std::cos(s*x) + std::sin(s*x)/x); + } + rhs[1](i,j,k) = rhs_val; } } @@ -251,7 +389,18 @@ void actual_init_prob_cyl1d (int i, int j, int k, } else { cce = std::sin(s*x)*s/x + std::cos(s*x)*s*s; } - rhs[2](i,j,k) = alpha*cce + beta*sol[2](i,j,k); + Real alpha_loc = alpha; + Real dalpha = Real(0.0); + if (alpha_node) { + alpha_loc = alpha_val(x); + dalpha = dalpha_val(x); + } + Real rhs_val = alpha_loc*cce + beta*sol[2](i,j,k); + if (alpha_node) { + // variable-alpha correction = +alpha'(r)*(curl E)_theta = -alpha'(r)*dE_z/dr + rhs_val += dalpha * s * std::sin(s*x); + } + rhs[2](i,j,k) = rhs_val; } } } diff --git a/Tests/LinearSolvers/CurlCurl/inputs b/Tests/LinearSolvers/CurlCurl/inputs index c3096d9a55..541b71c1c3 100644 --- a/Tests/LinearSolvers/CurlCurl/inputs +++ b/Tests/LinearSolvers/CurlCurl/inputs @@ -5,7 +5,7 @@ n_cell = 128 max_grid_size = 64 verbose = 2 -bottom_verbose = 2 +bottom_verbose = 0 alpha = 1.0 beta_factor = 0.01 diff --git a/Tests/LinearSolvers/CurlCurl/main.cpp b/Tests/LinearSolvers/CurlCurl/main.cpp index 76ca0b2e76..ac3e8a3e90 100644 --- a/Tests/LinearSolvers/CurlCurl/main.cpp +++ b/Tests/LinearSolvers/CurlCurl/main.cpp @@ -2,6 +2,9 @@ #include #include +#include + +using namespace amrex; int main (int argc, char* argv[]) { @@ -9,8 +12,36 @@ int main (int argc, char* argv[]) { BL_PROFILE("main"); - MyTest mytest; - mytest.solve(); + { + ParmParse pp; + Vector variable_alpha; + Vector variable_beta; + if (pp.contains("variable_alpha")) { + int t; + pp.get("variable_alpha", t); + variable_alpha.push_back(t); + } else { + variable_alpha.push_back(0); + variable_alpha.push_back(1); + } + if (pp.contains("variable_beta")) { + int t; + pp.get("variable_beta", t); + variable_beta.push_back(t); + } else { + variable_beta.push_back(0); + variable_beta.push_back(1); + } + for (auto alpha : variable_alpha) { + for (auto beta : variable_beta) { + pp.add("variable_alpha", alpha); + pp.add("variable_beta", beta); + amrex::Print() << "\nTesting variable_alpha = " << alpha << " variable_beta = " << beta << "\n"; + MyTest mytest; + mytest.solve(); + } + } + } } amrex::Finalize();