diff --git a/Src/AmrCore/AMReX_ErrorList.H b/Src/AmrCore/AMReX_ErrorList.H index 9bdbea056bd..3ca69b39a84 100644 --- a/Src/AmrCore/AMReX_ErrorList.H +++ b/Src/AmrCore/AMReX_ErrorList.H @@ -8,9 +8,13 @@ #include #include #include +#include +#include +#include namespace amrex { + extern "C" { @@ -372,6 +376,84 @@ private: std::ostream& operator << (std::ostream& os, const ErrorList& elst); + struct AMRErrorTagInfo + { + int m_max_level = 100000; + Real m_min_time = std::numeric_limits::lowest(); + Real m_max_time = std::numeric_limits::max(); + RealBox m_realbox; + + AMRErrorTagInfo& SetMaxLevel (int max_level) noexcept { + m_max_level = max_level; + return *this; + } + AMRErrorTagInfo& SetMinTime (amrex::Real min_time) noexcept { + m_min_time = min_time; + return *this; + } + AMRErrorTagInfo& SetMaxTime (amrex::Real max_time) noexcept { + m_max_time = max_time; + return *this; + } + AMRErrorTagInfo& SetRealBox (const amrex::RealBox& realbox) noexcept { + m_realbox = realbox; + return *this; + } + }; + + class AMRErrorTag + { + public: + + enum TEST {GRAD=0, LESS, GREATER, VORT, BOX, USER}; + + struct UserFunc + { + virtual void operator() (const amrex::Box& bx, + amrex::Array4 const& dat, + amrex::Array4 const& tag, + amrex::Real time, + int level, + char tagval, + char clearval) = 0; + }; + + explicit AMRErrorTag (const AMRErrorTagInfo& info = AMRErrorTagInfo()) noexcept + : m_test(BOX), m_field(std::string()), m_info(info) {m_ngrow = SetNGrow();} + + AMRErrorTag (amrex::Real value, + AMRErrorTag::TEST test, + const std::string& field, + const AMRErrorTagInfo& info = AMRErrorTagInfo()) noexcept + : m_value(value), m_test(test), m_field(field), m_info(info) {m_ngrow = SetNGrow();} + + AMRErrorTag (AMRErrorTag::UserFunc* userfunc, + const std::string& field, + int ngrow, + const AMRErrorTagInfo& info = AMRErrorTagInfo()) noexcept + : m_userfunc(userfunc), m_field(field), m_ngrow(ngrow), m_info(info) {} + + virtual void operator() (amrex::TagBoxArray& tb, + const amrex::MultiFab* mf, + char clearval, + char tagval, + amrex::Real time, + int level, + const amrex::Geometry& geom) const noexcept; + + int NGrow() const noexcept {return m_ngrow;} + const std::string& Field () const noexcept {return m_field;} + + protected: + int SetNGrow () const noexcept; + + Real m_value; + TEST m_test; + UserFunc* m_userfunc = nullptr; + std::string m_field; + AMRErrorTagInfo m_info; + int m_ngrow; + }; } #endif diff --git a/Src/AmrCore/AMReX_ErrorList.cpp b/Src/AmrCore/AMReX_ErrorList.cpp index a41bc1a199f..564aadf0a09 100644 --- a/Src/AmrCore/AMReX_ErrorList.cpp +++ b/Src/AmrCore/AMReX_ErrorList.cpp @@ -228,4 +228,190 @@ operator << (std::ostream& os, return os; } + static + void + AMRErrorTag_GRAD(const Box& bx, + Array4 const& dat, + Array4 const& tag, + Real threshold, + char tagval) + { + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + auto ax = amrex::Math::abs(dat(i+1,j,k) - dat(i,j,k)); + ax = amrex::max(ax,amrex::Math::abs(dat(i,j,k) - dat(i-1,j,k))); +#if AMREX_SPACEDIM == 1 + if (ax >= threshold) tag(i,j,k) = tagval; +#else + auto ay = amrex::Math::abs(dat(i,j+1,k) - dat(i,j,k)); + ay = amrex::max(ay,amrex::Math::abs(dat(i,j,k) - dat(i,j-1,k))); +#if AMREX_SPACEDIM > 2 + auto az = amrex::Math::abs(dat(i,j,k+1) - dat(i,j,k)); + az = amrex::max(az,amrex::Math::abs(dat(i,j,k) - dat(i,j,k-1))); +#endif + if (amrex::max(D_DECL(ax,ay,az)) >= threshold) { + tag(i,j,k) = tagval; + } +#endif + }); + } + + int + AMRErrorTag::SetNGrow () const noexcept + { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_test != USER, "Do not call SetNGrow with USER test"); + static std::map ng = { {GRAD,1}, {LESS,0}, {GREATER,0}, {VORT,0}, {BOX,0} }; + return ng[m_test]; + } + + static + void + AMRErrorTag_LESS(const Box& bx, + Array4 const& dat, + Array4 const& tag, + Real threshold, + char tagval) noexcept + { + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + if (dat(i,j,k) <= threshold) { + tag(i,j,k) = tagval; + } + }); + } + + static + void + AMRErrorTag_GREATER(const Box& bx, + Array4 const& dat, + Array4 const& tag, + Real threshold, + char tagval) noexcept + { + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + if (dat(i,j,k) >= threshold) { + tag(i,j,k) = tagval; + } + }); + } + + static + void + AMRErrorTag_BOX(const Box& bx, + Array4 const& tag, + const RealBox& tag_rb, + const Geometry& geom, + char tagval) noexcept + { + auto plo = geom.ProbLoArray(); + auto dx = geom.CellSizeArray(); + class RealBox trb(bx,dx.data(),plo.data()); + if (tag_rb.intersects(trb)) + { + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + GpuArray pt = {D_DECL(plo[0]+(i+0.5)*dx[0], + plo[1]+(j+0.5)*dx[1], + plo[2]+(k+0.5)*dx[2])}; + if (tag_rb.contains(pt.data())) { + tag(i,j,k) = tagval; + } + }); + } + } + + static + void + AMRErrorTag_VORT(const Box& bx, + Array4 const& dat, + Array4 const& tag, + int level, + Real threshold, + char tagval) noexcept + { + const Real fac = threshold * std::pow(2,level); + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + if (dat(i,j,k) >= fac) { + tag(i,j,k) = tagval; + } + }); + } + + void + AMRErrorTag::operator() (TagBoxArray& tba, + const MultiFab* mf, + char clearval, + char tagval, + Real time, + int level, + const Geometry& geom) const noexcept + { + BL_PROFILE("AMRErrorTag::operator()"); + + if (m_test == USER) + { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_userfunc!=nullptr,"UserFunc not properly set in AMRErrorTag"); + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(tba,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const auto& bx = mfi.tilebox(); + auto const& dat = mf->array(mfi); + auto tag = tba.array(mfi); + (*m_userfunc)(bx,dat,tag,time,level,tagval,clearval); + } + } + else + { + if ((level < m_info.m_max_level) && + (time >= m_info.m_min_time ) && + (time <= m_info.m_max_time ) ) + { + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(tba,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const auto& bx = mfi.tilebox(); + auto const& dat = mf->array(mfi); + auto tag = tba.array(mfi); + + if (m_test == BOX) + { + AMRErrorTag_BOX(bx, tag, m_info.m_realbox, geom, tagval); + } + else if (m_test == GRAD) + { + AMRErrorTag_GRAD(bx, dat, tag, m_value, tagval); + } + else if (m_test == LESS) + { + AMRErrorTag_LESS(bx, dat, tag, m_value, tagval); + } + else if (m_test == GREATER) + { + AMRErrorTag_GREATER(bx, dat, tag, m_value, tagval); + } + else if (m_test == VORT) + { + AMRErrorTag_VORT(bx, dat, tag, level, m_value, tagval); + } + else + { + Abort("Bad AMRErrorTag test flag"); + } + } + } + } + } }