Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
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
3 changes: 3 additions & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@ set(COMPONENT_NAME azplugins)
set(_${COMPONENT_NAME}_sources
ConstantFlow.cc
export_ImagePotentialBondHarmonic.cc
export_PerturbedLennardJonesEvap.cc
export_VariantInterpolated.cc
module.cc
ParabolicFlow.cc
VelocityCompute.cc
Expand Down Expand Up @@ -32,6 +34,7 @@ set(python_files
external.py
flow.py
pair.py
variant.py
wall.py
)

Expand Down
92 changes: 92 additions & 0 deletions src/LinearInterpolator1D.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,92 @@
// Copyright (c) 2018-2020, Michael P. Howard
// Copyright (c) 2021-2025, Auburn University
// Part of azplugins, released under the BSD 3-Clause License.

#ifndef AZPLUGINS_LINEAR_INTERPOLATOR_1D_H_
#define AZPLUGINS_LINEAR_INTERPOLATOR_1D_H_

#include <cassert>
#include <cmath>
#include <cstdint>

#include "hoomd/HOOMDMath.h"

#if defined(__HIPCC__) || defined(__CUDACC__)
#define AZPLUGINS_HOSTDEVICE __host__ __device__
#define AZPLUGINS_FORCEINLINE __forceinline__
#else
#define AZPLUGINS_HOSTDEVICE
#define AZPLUGINS_FORCEINLINE inline
#endif

namespace hoomd
{
namespace azplugins
{

template<typename T> class LinearInterpolator1D
{
public:
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE LinearInterpolator1D()
: m_data(nullptr), m_lo(Scalar(0)), m_hi(Scalar(0)), m_dx(Scalar(0)), m_n(0)
{
}

//! Constructor
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE
LinearInterpolator1D(const T* data, unsigned int n, Scalar lo, Scalar hi)
: m_data(data), m_lo(lo), m_hi(hi), m_n(n)
{
assert(n >= 2);
m_dx = (m_hi - m_lo) / Scalar(n - 1);
}

//! Implement piecewise linear interpolation
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar operator()(Scalar x) const
{
const Scalar f = (x - m_lo) / m_dx;
int b = static_cast<int>(std::floor(f));

// If exactly at the top boundary, shift into the last valid cell
if (f == Scalar(m_n - 1) && x == m_hi)
{
--b;
}

assert(b >= 0);
assert(b < static_cast<int>(m_n - 1));
const Scalar frac = f - Scalar(b);

const Scalar c0 = m_data[b];
const Scalar c1 = m_data[b + 1];

// Linear interpolation
return c0 * (Scalar(1) - frac) + c1 * frac;
}

//! Lower bound
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getLo() const
{
return m_lo;
}
//! Upper bound
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getHi() const
{
return m_hi;
}

private:
const T* m_data;
Scalar m_lo;
Scalar m_hi;
Scalar m_dx;
unsigned int m_n;
};

} // namespace azplugins
} // namespace hoomd

#undef AZPLUGINS_HOSTDEVICE
#undef AZPLUGINS_FORCEINLINE

#endif // AZPLUGINS_LINEAR_INTERPOLATOR_1D_H_
130 changes: 130 additions & 0 deletions src/LinearInterpolator2D.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,130 @@
// Copyright (c) 2018-2020, Michael P. Howard
// Copyright (c) 2021-2025, Auburn University
// Part of azplugins, released under the BSD 3-Clause License.

#ifndef AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_
#define AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_

#include <cassert>
#include <cmath>
#include <cstdint>

#include "hoomd/HOOMDMath.h"
#include "hoomd/Index1D.h"

#if defined(__HIPCC__) || defined(__CUDACC__)
#define AZPLUGINS_HOSTDEVICE __host__ __device__
#define AZPLUGINS_FORCEINLINE __forceinline__
#else
#define AZPLUGINS_HOSTDEVICE
#define AZPLUGINS_FORCEINLINE inline
#endif

namespace hoomd
{
namespace azplugins
{

template<typename T> class LinearInterpolator2D
{
public:
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE LinearInterpolator2D() : m_data(nullptr), m_indexer()
{
m_n[0] = 0;
m_n[1] = 0;
for (int d = 0; d < 2; ++d)
{
m_lo[d] = Scalar(0);
m_hi[d] = Scalar(0);
m_dx[d] = Scalar(0);
}
}

//! Constructor
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE
LinearInterpolator2D(const T* data, const unsigned int* n, const Scalar* lo, const Scalar* hi)
: m_data(data), m_indexer(n[0], n[1])
{
for (int d = 0; d < 2; ++d)
{
const unsigned int nd = n[d];
assert(nd >= 2);

m_n[d] = nd;
m_lo[d] = lo[d];
m_hi[d] = hi[d];
m_dx[d] = (m_hi[d] - m_lo[d]) / Scalar(nd - 1);
}
}

//! Interpolate at (x0, x1)
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar operator()(Scalar x0, Scalar x1) const
{
const Scalar x[2] = {x0, x1};

int bin[2];
Scalar frac[2];

for (int d = 0; d < 2; ++d)
{
const unsigned int nd = m_n[d];
const Scalar f = (x[d] - m_lo[d]) / m_dx[d];
int b = static_cast<int>(std::floor(f));

// If exactly at the top boundary, shift into the last valid cell
if (f == Scalar(nd - 1) && x[d] == m_hi[d])
{
--b;
}

assert(b >= 0);
assert(b < static_cast<int>(nd - 1));
bin[d] = b;
frac[d] = f - Scalar(b);
}

const unsigned int b0 = static_cast<unsigned int>(bin[0]);
const unsigned int b1 = static_cast<unsigned int>(bin[1]);
const Scalar c00 = m_data[m_indexer(b0, b1)];
const Scalar c01 = m_data[m_indexer(b0, b1 + 1)];
const Scalar c10 = m_data[m_indexer(b0 + 1, b1)];
const Scalar c11 = m_data[m_indexer(b0 + 1, b1 + 1)];

// Bilinear interpolation
const Scalar t0 = frac[0];
const Scalar omt0 = Scalar(1) - t0;
const Scalar c0 = c00 * omt0 + c10 * t0;
const Scalar c1 = c01 * omt0 + c11 * t0;

const Scalar t1 = frac[1];
return c0 * (Scalar(1) - t1) + c1 * t1;
}

//! Lower bound
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getLo(int dim) const
{
return m_lo[dim];
}

//! Upper bound
AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getHi(int dim) const
{
return m_hi[dim];
}

private:
const T* m_data;
Scalar m_lo[2];
Scalar m_hi[2];
Scalar m_dx[2];
unsigned int m_n[2];
Index2D m_indexer;
};

} // namespace azplugins
} // namespace hoomd

#undef AZPLUGINS_HOSTDEVICE
#undef AZPLUGINS_FORCEINLINE

#endif // AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_
Loading
Loading