Skip to content

Commit

Permalink
begin 2dof test case
Browse files Browse the repository at this point in the history 
- fixed tensor slicing
- backend blas abstraction
- begin newmark method implementation
  • Loading branch information
@samayala22
samayala22 committed Sep 19, 2024
1 parent 02b4da6 commit 0adcf51
Show file tree
Hide file tree
Showing 8 changed files with 344 additions and 94 deletions.
185 changes: 96 additions & 89 deletions headeronly/tinyview.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ namespace tiny {

using Range = std::array<int, 2>;

constexpr Range all{0,-1};
constexpr Range all{0, -1};

template<typename... Args>
struct CountRanges;
Expand All @@ -31,115 +31,122 @@ struct CountRanges<First, Rest...> {

template<typename T, std::size_t N>
class View {
public:
constexpr View(T* const ptr, const std::array<std::size_t, N>& dims) : _ptr(ptr), _dims(dims), _strides(compute_strides(dims)) {}
constexpr View(T* const ptr, const std::array<std::size_t, N>& dims, const std::array<std::size_t, N>& strides) : _ptr(ptr), _dims(dims), _strides(strides) {}
constexpr View(const View& other) : _ptr(other._ptr), _dims(other._dims), _strides(other._strides) {}
template<typename... Idx>
TINY_INLINE constexpr T& operator()(std::size_t first, Idx... idx) {return _ptr[compute_index(first, idx...)];}

template<typename... Idx>
TINY_INLINE constexpr T* ptr(std::size_t first, Idx... idx) {return _ptr + compute_index(first, idx...);}

template<typename... Args>
TINY_INLINE constexpr auto view(Args... args) {
constexpr std::size_t M = CountRanges<Args...>::value;
static_assert(sizeof...(args) == N, "The number of indices must match the dimension N.");
static_assert(M <= N, "Too many ranges provided compared to the view's dimensionality");

T* newPtr = _ptr;
std::array<std::size_t, M> newDims, newStrides;
std::size_t newDimIndex = 0;

size_t argIndex = 0;
([&](auto& arg) {
if constexpr (std::is_same<std::decay_t<decltype(arg)>, Range>::value) {
constexpr size_t first = arg[0];
constexpr size_t last = (arg[1] < 0) ? _dims[argIndex] + arg[1] + 1 : arg[1];
assert((first >= 0) && (first < _dims[argIndex]));
assert((last >= 0) && (last < _dims[argIndex]));
assert(last - first > 0);
newPtr += first * _strides[argIndex];
newDims[newDimIndex] = last - first;
newStrides[newDimIndex] = _strides[argIndex];
newDimIndex++;
} else if constexpr (std::is_integral<std::decay_t<decltype(arg)>>::value) {
constexpr size_t real_arg = (arg < 0) ? _dims[argIndex] + arg + 1 : arg;
newPtr += real_arg * _strides[argIndex];
}
argIndex++;
}(args), ...);

return View<T, M>(newPtr, newDims, newStrides);
}

TINY_INLINE constexpr T* data() {return _ptr;}
TINY_INLINE constexpr View& operator=(const View& other) = default;

TINY_INLINE constexpr const std::size_t& dims(std::size_t idx) const {return _dims[idx];}
TINY_INLINE constexpr const std::size_t& strides(std::size_t idx) const {return _strides[idx];}

private:
template<typename... Idx>
TINY_INLINE constexpr std::size_t compute_index(std::size_t first, Idx... idx) {
static_assert(sizeof...(idx) == N-1, "The number of indices must match the dimension N.");
static_assert((std::is_integral<std::decay_t<Idx>>::value && ...), "All indices must be integers");

#ifndef NDEBUG
std::size_t ii = 1;
assert(((idx < _dims[ii++]) && ...));
#endif

std::size_t index = first;
std::size_t i = 1;
((index += idx * _strides[i++]), ...);
return index;
}

TINY_INLINE constexpr static std::array<std::size_t, N> compute_strides(const std::array<std::size_t, N>& dims) {
std::array<std::size_t, N> s{};
s[0] = 1;
for (std::size_t i = 1; i < N; ++i) {
s[i] = s[i - 1] * dims[i - 1];
public:
constexpr View(T* const ptr, const std::array<std::size_t, N>& dims)
: _ptr(ptr), _dims(dims), _strides(compute_strides(dims)) {}

constexpr View(T* const ptr, const std::array<std::size_t, N>& dims, const std::array<std::size_t, N>& strides)
: _ptr(ptr), _dims(dims), _strides(strides) {}

constexpr View(const View& other)
: _ptr(other._ptr), _dims(other._dims), _strides(other._strides) {}

template<typename... Idx>
TINY_INLINE T& operator()(std::size_t first, Idx... idx) const { return _ptr[compute_index(first, idx...)]; }

template<typename... Idx>
TINY_INLINE T* ptr(std::size_t first, Idx... idx) const { return _ptr + compute_index(first, idx...); }

template<typename... Args>
TINY_INLINE auto view(Args... args) { // Removed 'constexpr'
constexpr std::size_t M = CountRanges<Args...>::value;
static_assert(sizeof...(args) == N, "The number of indices must match the dimension N.");
static_assert(M <= N, "Too many ranges provided compared to the view's dimensionality");

T* newPtr = _ptr;
std::array<std::size_t, M> newDims{};
std::array<std::size_t, M> newStrides{};
std::size_t newDimIndex = 0;

std::size_t argIndex = 0;
([&](auto& arg) {
if constexpr (std::is_same<std::decay_t<decltype(arg)>, Range>::value) {
std::size_t first = arg[0]; // Removed 'constexpr'
std::size_t last = (arg[1] < 0) ? _dims[argIndex] + arg[1] + 1 : arg[1];
assert((first >= 0) && (first < _dims[argIndex]));
assert((last >= 0) && (last <= _dims[argIndex])); // Changed '<' to '<=' for upper bound
assert(last - first > 0);
newPtr += first * _strides[argIndex];
newDims[newDimIndex] = last - first;
newStrides[newDimIndex] = _strides[argIndex];
newDimIndex++;
} else if constexpr (std::is_integral<std::decay_t<decltype(arg)>>::value) {
std::size_t real_arg = (arg < 0) ? _dims[argIndex] + arg : arg; // Removed '+1' as indexing starts from 0
assert((real_arg >= 0) && (real_arg < _dims[argIndex]));
newPtr += real_arg * _strides[argIndex];
}
return s;
argIndex++;
}(args), ...);

return View<T, M>(newPtr, newDims, newStrides);
}

TINY_INLINE T* data() const { return _ptr; }
TINY_INLINE View& operator=(const View& other) = default;

TINY_INLINE const std::size_t& dims(std::size_t idx) const { return _dims[idx]; }
TINY_INLINE const std::size_t& strides(std::size_t idx) const { return _strides[idx]; }

private:
template<typename... Idx>
std::size_t compute_index(std::size_t first, Idx... idx) const { // Removed 'constexpr'
static_assert(sizeof...(idx) == N-1, "The number of indices must match the dimension N.");
static_assert((std::is_integral<std::decay_t<Idx>>::value && ...), "All indices must be integers");

#ifndef NDEBUG
std::size_t ii = 0; // Changed from 1 to 0 to match zero-based indexing
((assert(idx >= 0 && idx < _dims[ii++])), ...);
#endif

std::size_t index = first;
std::size_t i = 1;
((index += idx * _strides[i++]), ...);
return index;
}

static std::array<std::size_t, N> compute_strides(const std::array<std::size_t, N>& dims) {
std::array<std::size_t, N> s{};
s[0] = 1;
for (std::size_t i = 1; i < N; ++i) {
s[i] = s[i - 1] * dims[i - 1];
}
return s;
}

T* const _ptr;
const std::array<std::size_t, N> _dims;
const std::array<std::size_t, N> _strides;
T* const _ptr;
const std::array<std::size_t, N> _dims;
const std::array<std::size_t, N> _strides;
};
} // namespace tiny

} // namespace tiny

// #include <cstdio>
// using namespace tiny;

// int main() {
// std::size_t n = 3;
// float a[27] = {};
// constexpr int n = 3;
// float a[n * n * n] = {};

// View<float,3> av{&a[0], {n,n,n}};

// foo2(av);
// for (int i = 0; i < n * n * n; i++) {
// a[i] = static_cast<float>(i);
// }

// View<float, 2> bv = av.view(all, -2, Range{0,3});
// View<float, 3> av{&a[0], {static_cast<std::size_t>(n), static_cast<std::size_t>(n), static_cast<std::size_t>(n)}};

// View<float, 2> bv = av.view(all, -2, Range{0, 3});

// for (int j = 0; j < bv.dims(1); j++) {
// for (int i = 0; i < bv.dims(0); i++) {
// std::printf("%f ", bv(i,j));
// for (std::size_t j = 0; j < bv.dims(1); j++) {
// for (std::size_t i = 0; i < bv.dims(0); i++) {
// std::printf("%f ", bv(i, j));
// }
// std::printf("\n");
// }

// assert(bv() - av() == 3);
// assert(bv.dims(0) == 3);
// assert(bv.dims(1) == 3);
// assert(bv(0,0) == 3.f);
// assert(bv(2,1) == 14.f);
// assert(bv(2,2) == 23.f);

// // for (int i = 0; i < 27; i++) {
// // std::printf("%f ", a[i]);
// // }
// return 0;
// }
107 changes: 107 additions & 0 deletions tests/uvlm_2dof.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
#include "vlm.hpp"
#include "vlm_utils.hpp"
#include "vlm_kinematics.hpp"

#include "tinycombination.hpp"

using namespace vlm;

class NewmarkBeta {
public:
NewmarkBeta(std::unique_ptr<Memory> memory, const f32 beta = 0.25f, const f32 gamma = 0.5f) : m_memory(std::move(memory)), m_beta(beta), m_gamma(gamma) {};
~NewmarkBeta() = default;

void run(View<f32, Tensor<2>>& M, View<f32, Tensor<2>>& C, View<f32, Tensor<2>>& K, View<f32, Tensor<1>>& F, View<f32, Tensor<1>>& u0, View<f32, Tensor<1>>& v0, View<f32, Tensor<1>>& t);

private:
std::unique_ptr<BLAS> m_blas;
std::unique_ptr<Memory> m_memory;
const f32 m_beta;
const f32 m_gamma;

Buffer<f32, MemoryLocation::Device, Tensor<2>> K_eff{*m_memory}; // effective stiffness
Buffer<f32, MemoryLocation::Device, Tensor<1>> a0{*m_memory}; // initial acceleration
Buffer<f32, MemoryLocation::Device, Tensor<1>> du{*m_memory}; // incremental displacement
Buffer<f32, MemoryLocation::Device, Tensor<1>> factor{*m_memory}; // intermediary vector
Buffer<f32, MemoryLocation::Device, Tensor<2>> u{*m_memory}; // dof x tsteps position history
Buffer<f32, MemoryLocation::Device, Tensor<2>> v{*m_memory}; // dof x tsteps velocity history
Buffer<f32, MemoryLocation::Device, Tensor<2>> a{*m_memory}; // dof x tsteps acceleration history
};

void NewmarkBeta::run(View<f32, Tensor<2>>& M, View<f32, Tensor<2>>& C, View<f32, Tensor<2>>& K, View<f32, Tensor<1>>& F, View<f32, Tensor<1>>& u0, View<f32, Tensor<1>>& v0, View<f32, Tensor<1>>& t) {
assert(M.layout.shape(0) == C.layout.shape(0));
assert(M.layout.shape(1) == C.layout.shape(1));
assert(C.layout.shape(0) == K.layout.shape(0));
assert(C.layout.shape(1) == K.layout.shape(1));
assert(K.layout.shape(0) == F.layout.shape(0));
assert(u0.layout.shape(0) == F.layout.shape(0));
assert(v0.layout.shape(0) == F.layout.shape(0));

Tensor<2> time_series{{F.layout.shape(0), t.layout.shape(0)}}; // dofs x timesteps

// TODO: preallocate a number of timesteps and only resize when necessary
u.dealloc();
v.dealloc();
a.dealloc();

u.alloc(time_series);
v.alloc(time_series);
a.alloc(time_series);

m_memory->copy(MemoryTransfer::DeviceToDevice, u.d_view().ptr, u0.ptr, u0.size_bytes());
m_memory->copy(MemoryTransfer::DeviceToDevice, v.d_view().ptr, v0.ptr, v0.size_bytes());

// a[:,0] = np.linalg.solve(M, F[0] - C @ v0 - K @ x0)
View<f32, Tensor<1>> a0_col0 = a.d_view().layout.slice(a.d_view().ptr, all, 0);
m_memory->copy(MemoryTransfer::DeviceToDevice, a0_col0.ptr, F.ptr, F.size_bytes());
m_blas->gemv(-1.0f, C, v0, 1.0f, a0_col0);
m_blas->gemv(-1.0f, K, u0, 1.0f, a0_col0);
}

int main() {
const u64 ni = 20;
const u64 nj = 5;
// vlm::Executor::instance(1);
const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_" + std::to_string(ni) + "x" + std::to_string(nj) + ".x"};
const std::vector<std::string> backends = {"cpu"};

auto solvers = tiny::make_combination(meshes, backends);

// Geometry
const f32 b = 0.5f; // half chord

// Define simulation length
const f32 cycles = 10.0f;
const f32 u_inf = 1.0f; // freestream velocity
const f32 amplitude = 0.1f; // amplitude of the wing motion
const f32 k = 0.5; // reduced frequency
const f32 omega = k * 2.0f * u_inf / (2*b);
const f32 t_final = cycles * 2.0f * PI_f / omega; // 4 periods
//const f32 t_final = 5.0f;

Kinematics kinematics{};

const f32 initial_angle = 0.0f;

const auto initial_pose = rotation_matrix(
linalg::alias::float3{0.0f, 0.0f, 0.0f}, // take into account quarter chord panel offset
linalg::alias::float3{0.0f, 1.0f, 0.0f},
to_radians(initial_angle)
);

// Sudden acceleration
const f32 alpha = to_radians(5.0f);
kinematics.add([=](const fwd::Float& t) {
return translation_matrix<fwd::Float>({
-u_inf*std::cos(alpha)*t,
0.0f,
-u_inf*std::sin(alpha)*t
});
});

for (const auto& [mesh_name, backend_name] : solvers) {
UVLM simulation{backend_name, {mesh_name}};
simulation.run({kinematics}, {initial_pose}, t_final);
}
return 0;
}
9 changes: 9 additions & 0 deletions vlm/backends/cpu/include/vlm_backend_cpu.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,4 +35,13 @@ class BackendCPU final : public Backend {
f32 mesh_area(const View<f32, SingleSurface>& areas) override;
};

class BLAS_CPU final : public BLAS {
public:
BLAS_CPU() = default;
~BLAS_CPU() = default;

void gemv(const f32 alpha, const View<f32, Tensor<2>>& A, const View<f32, Tensor<1>>& x, const f32 beta, View<f32, Tensor<1>>& y, Trans trans = Trans::No) override;
void gemm(const f32 alpha, const View<f32, Tensor<2>>& A, const View<f32, Tensor<2>>& B, const f32 beta, View<f32, Tensor<2>>& C, Trans trans_a = Trans::No, Trans trans_b = Trans::No) override;
};

} // namespace vlm
Loading

0 comments on commit 0adcf51

Please sign in to comment.