Skip to content

Commit

Permalink
Use linear interpolation in DragForcing (#1323)
Browse files Browse the repository at this point in the history
  • Loading branch information
@marchdf
marchdf authored Nov 4, 2024
1 parent cd0c0f1 commit 2efa3f8
Show file tree
Hide file tree
Showing 3 changed files with 148 additions and 32 deletions.
33 changes: 16 additions & 17 deletions amr-wind/equation_systems/icns/source_terms/DragForcing.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include "AMReX_Gpu.H"
#include "AMReX_Random.H"
#include "amr-wind/wind_energy/ABL.H"
#include "amr-wind/utilities/linear_interpolation.H"

namespace amr_wind::pde::icns {

Expand Down Expand Up @@ -115,23 +116,21 @@ void DragForcing::operator()(
yi_end = sponge_north * std::max(yi_end, 0.0);
ystart_damping = sponge_strength * yi_start * yi_start;
yend_damping = sponge_strength * yi_end * yi_end;
amrex::Real spongeVelX = 0.0;
amrex::Real spongeVelY = 0.0;
amrex::Real spongeVelZ = 0.0;
amrex::Real residual = 1000;
amrex::Real height_error = 0.0;
for (unsigned ii = 0; ii < vsize; ++ii) {
height_error = std::abs(z - device_vel_ht[ii]);
if (height_error < residual) {
residual = height_error;
const unsigned ix = 3 * ii;
const unsigned iy = 3 * ii + 1;
const unsigned iz = 3 * ii + 2;
spongeVelX = device_vel_vals[ix];
spongeVelY = device_vel_vals[iy];
spongeVelZ = device_vel_vals[iz];
}
}

const auto idx =
interp::bisection_search(device_vel_ht, device_vel_ht + vsize, z);
const amrex::Real spongeVelX =
(vsize > 0) ? interp::linear_impl(
device_vel_ht, device_vel_vals, z, idx, 3, 0)
: 0.0;
const amrex::Real spongeVelY =
(vsize > 0) ? interp::linear_impl(
device_vel_ht, device_vel_vals, z, idx, 3, 1)
: 0.0;
const amrex::Real spongeVelZ =
(vsize > 0) ? interp::linear_impl(
device_vel_ht, device_vel_vals, z, idx, 3, 2)
: 0.0;
amrex::Real Dxz = 0.0;
amrex::Real Dyz = 0.0;
amrex::Real bc_forcing_x = 0;
Expand Down
82 changes: 67 additions & 15 deletions amr-wind/utilities/linear_interpolation.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,23 @@ bisection_search(const It begin, const It end, const T& x)
return idx;
}

template <typename It, typename T>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Index
nearest_search(const It begin, const It end, const T& x)
{

auto idx = bisection_search(begin, end, x);
if ((idx.lim == Limits::LOWLIM) || (idx.lim == Limits::UPLIM)) {
return idx;
}

if ((x - begin[idx.idx]) > (begin[idx.idx + 1] - x)) {
idx.idx += 1;
}

return idx;
}

template <typename It, typename T>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Index
find_index(const It begin, const It end, const T& x, const int hint = 1)
Expand All @@ -82,36 +99,63 @@ find_index(const It begin, const It end, const T& x, const int hint = 1)
template <typename C1, typename C2>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
typename std::iterator_traits<C2>::value_type
linear(
linear_impl(
const C1 xbegin,
const C1 xend,
const C2 yinp,
const typename std::iterator_traits<C1>::value_type& xout)
const typename std::iterator_traits<C1>::value_type& xout,
const Index& idx,
const int ncomp = 1,
const int comp = 0)
{
using DType1 = typename std::iterator_traits<C1>::value_type;
const auto idx = bisection_search(xbegin, xend, xout);

if ((idx.lim == Limits::LOWLIM) || (idx.lim == Limits::UPLIM)) {
return yinp[idx.idx];
return yinp[ncomp * idx.idx + comp];
}
static constexpr DType1 eps = 1.0e-8;
const int j = idx.idx;
const auto denom = (xbegin[j + 1] - xbegin[j]);
const auto facR = (denom > eps) ? ((xout - xbegin[j]) / denom) : 1.0;
const auto facL = static_cast<DType1>(1.0) - facR;
return facL * yinp[j] + facR * yinp[j + 1];
return facL * yinp[ncomp * j + comp] + facR * yinp[ncomp * (j + 1) + comp];
}

template <typename C1, typename C2>
inline typename C2::value_type
linear(const C1& xinp, const C2& yinp, const typename C1::value_type& xout)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
typename std::iterator_traits<C2>::value_type
linear(
const C1 xbegin,
const C1 xend,
const C2 yinp,
const typename std::iterator_traits<C1>::value_type& xout,
const int ncomp = 1,
const int comp = 0)
{
const auto idx = bisection_search(xbegin, xend, xout);
return linear_impl(xbegin, yinp, xout, idx, ncomp, comp);
}

template <typename C1, typename C2>
inline typename C2::value_type linear(
const C1& xinp,
const C2& yinp,
const typename C1::value_type& xout,
const int ncomp = 1,
const int comp = 0)
{
return linear(xinp.data(), (xinp.data() + xinp.size()), yinp.data(), xout);
return linear(
xinp.data(), (xinp.data() + xinp.size()), yinp.data(), xout, ncomp,
comp);
}

template <typename C1, typename C2>
inline void
linear_monotonic(const C1& xinp, const C2& yinp, const C1& xout, C2& yout)
inline void linear_monotonic(
const C1& xinp,
const C2& yinp,
const C1& xout,
C2& yout,
const int ncomp = 1,
const int comp = 0)
{
static constexpr typename C1::value_type eps = 1.0e-8;
AMREX_ASSERT(xinp.size() == yinp.size());
Expand All @@ -125,20 +169,28 @@ linear_monotonic(const C1& xinp, const C2& yinp, const C1& xout, C2& yout)
find_index(xinp.data(), (xinp.data() + xinp.size()), x, hint);

if ((idx.lim == Limits::LOWLIM) || (idx.lim == Limits::UPLIM)) {
yout[i] = yinp[idx.idx];
yout[i] = yinp[ncomp * idx.idx + comp];
} else if (idx.lim == Limits::VALID) {
int j = idx.idx;
const auto denom = (xinp[j + 1] - xinp[j]);
const auto facR = (denom > eps) ? ((x - xinp[j]) / denom) : 1.0;
const auto facL = static_cast<typename C1::value_type>(1.0) - facR;
yout[i] = facL * yinp[j] + facR * yinp[j + 1];
yout[i] = facL * yinp[ncomp * j + comp] +
facR * yinp[ncomp * (j + 1) + comp];
;
}
hint = idx.idx + 1;
}
}

template <typename C1, typename C2>
inline void linear(const C1& xinp, const C2& yinp, const C1& xout, C2& yout)
inline void linear(
const C1& xinp,
const C2& yinp,
const C1& xout,
C2& yout,
const int ncomp = 1,
const int comp = 0)
{
AMREX_ASSERT(
static_cast<amrex::Long>(xinp.size()) ==
Expand All @@ -149,7 +201,7 @@ inline void linear(const C1& xinp, const C2& yinp, const C1& xout, C2& yout)

int npts = xout.size();
for (int i = 0; i < npts; ++i) {
yout[i] = linear(xinp, yinp, xout[i]);
yout[i] = linear(xinp, yinp, xout[i], ncomp, comp);
}
}

Expand Down
65 changes: 65 additions & 0 deletions unit_tests/utilities/test_linear_interpolation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,46 @@ TEST(LinearInterpolation, bisection_search)
}
}

TEST(LinearInterpolation, nearest_search)
{
namespace interp = amr_wind::interp;
std::vector<amrex::Real> xvec(10);
std::iota(xvec.begin(), xvec.end(), 0.0);

const auto* start = xvec.data();
const auto* end = xvec.data() + xvec.size();
{
const auto idx = interp::nearest_search(start, end, 5.0);
EXPECT_EQ(idx.idx, 5);
EXPECT_EQ(idx.lim, interp::Limits::VALID);
}
{
const auto idx = interp::nearest_search(start, end, 1.1);
EXPECT_EQ(idx.idx, 1);
EXPECT_EQ(idx.lim, interp::Limits::VALID);
}
{
const auto idx = interp::nearest_search(start, end, 1.6);
EXPECT_EQ(idx.idx, 2);
EXPECT_EQ(idx.lim, interp::Limits::VALID);
}
{
const auto idx = interp::nearest_search(start, end, 3.5);
EXPECT_EQ(idx.idx, 3);
EXPECT_EQ(idx.lim, interp::Limits::VALID);
}
{
const auto idx = interp::nearest_search(start, end, -1.0);
EXPECT_EQ(idx.idx, 0);
EXPECT_EQ(idx.lim, interp::Limits::LOWLIM);
}
{
const auto idx = interp::nearest_search(start, end, 9.1);
EXPECT_EQ(idx.idx, 9);
EXPECT_EQ(idx.lim, interp::Limits::UPLIM);
}
}

TEST(LinearInterpolation, find_index)
{
namespace interp = amr_wind::interp;
Expand Down Expand Up @@ -108,6 +148,31 @@ TEST(LinearInterpolation, lin_interp_single)
}
}

TEST(LinearInterpolation, lin_interp_single_multicomponent)
{
namespace interp = amr_wind::interp;

const int ncomp = 3;
std::vector<amrex::Real> xvec(10), yvec(ncomp * 10);
std::iota(xvec.begin(), xvec.end(), 0.0);
const amrex::Vector<amrex::Real> mult_facs = {
2.0 + 10.0 * amrex::Random(), 2.0 + 10.0 * amrex::Random(),
2.0 + 10.0 * amrex::Random()};
for (int i = 0; i < static_cast<int>(xvec.size()); i++) {
for (int n = 0; n < ncomp; n++) {
yvec[ncomp * i + n] = mult_facs[n] * xvec[i];
}
}

std::vector<amrex::Real> xtest{2.5, 4.5, 6.3, 8.8};
for (const auto& x : xtest) {
for (int n = 0; n < ncomp; n++) {
const auto y = interp::linear(xvec, yvec, x, 3, n);
EXPECT_NEAR(y, mult_facs[n] * x, 1.0e-12);
}
}
}

TEST(LinearInterpolation, lin_interp)
{
namespace interp = amr_wind::interp;
Expand Down

0 comments on commit 2efa3f8

Please sign in to comment.