Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add alternative error tagging class #1166

Merged
merged 2 commits into from
Jul 23, 2020
Merged

Add alternative error tagging class #1166

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
63d75e2
Add alternative error tagging class
drummerdoc Jul 17, 2020
c72ad8e
Address PR comments
drummerdoc Jul 23, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 82 additions & 0 deletions Src/AmrCore/AMReX_ErrorList.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,13 @@
#include <AMReX_Vector.H>
#include <AMReX_REAL.H>
#include <AMReX_ArrayLim.H>
#include <AMReX_MultiFab.H>
#include <AMReX_TagBox.H>
#include <AMReX_Geometry.H>

namespace amrex {


extern "C"
{

Expand Down Expand Up @@ -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<Real>::lowest();
Real m_max_time = std::numeric_limits<Real>::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 amrex::Real> const& dat,
amrex::Array4<char> 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
186 changes: 186 additions & 0 deletions Src/AmrCore/AMReX_ErrorList.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -228,4 +228,190 @@ operator << (std::ostream& os,
return os;
}

static
void
AMRErrorTag_GRAD(const Box& bx,
Array4<const Real> const& dat,
Array4<char> 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<TEST,int> ng = { {GRAD,1}, {LESS,0}, {GREATER,0}, {VORT,0}, {BOX,0} };
return ng[m_test];
}

static
void
AMRErrorTag_LESS(const Box& bx,
Array4<const Real> const& dat,
Array4<char> 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 Real> const& dat,
Array4<char> 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<char> 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<Real,AMREX_SPACEDIM> pt = {D_DECL(plo[0]+(i+0.5)*dx[0],
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@drummerdoc @WeiqunZhang Getting a AppleClang (SP) and DPC++ (DP) compile error on this of the form:

Src/AmrCore/AMReX_ErrorList.cpp:318:52: error: non-constant-expression cannot be narrowed from type 'double' to 'float' in initializer list [-Wc++11-narrowing]

in WarpX.

Copy link
Member

@ax3l ax3l Jul 23, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe we can cast the i, j, k with Real(i) etc. and add _rt in + 0.5 as + 0.5_rt to fix it.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's try this: #1180

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 Real> const& dat,
Array4<char> 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");
}
}
}
}
}
}