Skip to content

Commit

Permalink
slow progress
Browse files Browse the repository at this point in the history
  • Loading branch information
@samayala22
samayala22 committed Mar 27, 2024
1 parent ab06a00 commit a9d8923
Show file tree
Hide file tree
Showing 5 changed files with 35 additions and 43 deletions.
27 changes: 5 additions & 22 deletions headeronly/linalg.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -246,23 +246,6 @@ namespace linalg
};

// Small, fixed-size matrix type, consisting of exactly M rows and N columns of type T, stored in column-major order.
template<class T, int M> struct mat<T,M,1>
{
typedef vec<T,M> V;
V x;
constexpr mat() : x() {}
constexpr mat(const V & x_) : x(x_) {}
constexpr explicit mat(const T & s) : x(s) {}
constexpr explicit mat(const T * p) : x(p+M*0) {}
template<class U>
constexpr explicit mat(const mat<U,M,1> & m) : mat(V(m.x)) {}
constexpr vec<T,1> row(int i) const { return {x[i]}; }
constexpr const V & operator[] (int j) const { return x; }
constexpr V & operator[] (int j) { return x; }

template<class U, class=detail::conv_t<mat,U>> constexpr mat(const U & u) : mat(converter<mat,U>{}(u)) {}
template<class U, class=detail::conv_t<U,mat>> constexpr operator U () const { return converter<U,mat>{}(*this); }
};
template<class T, int M> struct mat<T,M,2>
{
typedef vec<T,M> V;
Expand All @@ -281,7 +264,7 @@ namespace linalg
template<class U, class=detail::conv_t<U,mat>> constexpr operator U () const { return converter<U,mat>{}(*this); }
};
template<class T, int M> struct mat<T,M,3>
{
{
typedef vec<T,M> V;
V x,y,z;
constexpr mat() : x(), y(), z() {}
Expand All @@ -292,8 +275,8 @@ namespace linalg
template<class U>
constexpr explicit mat(const mat<U,M,3> & m) : mat(V(m.x), V(m.y), V(m.z)) {}
constexpr vec<T,3> row(int i) const { return {x[i], y[i], z[i]}; }
constexpr const V & operator[] (int j) const { return j==0?x:j==1?y:z; }
constexpr V & operator[] (int j) { return j==0?x:j==1?y:z; }
constexpr const V & col (int j) const { return j==0?x:j==1?y:z; }
constexpr V & col (int j) { return j==0?x:j==1?y:z; }

template<class U, class=detail::conv_t<mat,U>> constexpr mat(const U & u) : mat(converter<mat,U>{}(u)) {}
template<class U, class=detail::conv_t<U,mat>> constexpr operator U () const { return converter<U,mat>{}(*this); }
Expand All @@ -310,8 +293,8 @@ namespace linalg
template<class U> // typecast matrix ?
constexpr explicit mat(const mat<U,M,4> & m) : mat(V(m.x), V(m.y), V(m.z), V(m.w)) {}
constexpr vec<T,4> row(int i) const { return {x[i], y[i], z[i], w[i]}; }
constexpr const V & operator[] (int j) const { return j==0?x:j==1?y:j==2?z:w; }
constexpr V & operator[] (int j) { return j==0?x:j==1?y:j==2?z:w; }
constexpr const V & col (int j) const { return j==0?x:j==1?y:j==2?z:w; }
constexpr V & col (int j) { return j==0?x:j==1?y:j==2?z:w; }

template<class U, class=detail::conv_t<mat,U>> constexpr mat(const U & u) : mat(converter<mat,U>{}(u)) {}
template<class U, class=detail::conv_t<U,mat>> constexpr operator U () const { return converter<U,mat>{}(*this); }
Expand Down
21 changes: 10 additions & 11 deletions tests/test_uvlm_theodorsen.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ int main() {
});
};

// For the moment take AoA=0 and U_inf=1
// For the moment take AoA=0 and U_inf=1
auto displacement_freestream = [&](f32 t) -> linalg::alias::float4x4 {
return linalg::translation_matrix(linalg::alias::float3{
-1.0f*t,
Expand All @@ -43,6 +43,11 @@ int main() {
return linalg::mul(displacement_freestream(t), displacement_wing(t));
};

auto relative_displacement = [&](f32 t0, f32 t1) -> linalg::alias::float4x4 {
return linalg::mul(displacement(t1), linalg::inverse(displacement(t0)));
};

// Magnitude of the velocity of a vertex at time t calculated by finite difference
auto absolute_velocity = [&](f32 t, const linalg::alias::float4& vertex) -> f32 {
return linalg::length((linalg::mul(displacement(t+EPS_sqrt_f), vertex)-linalg::mul(displacement(t), vertex))/EPS_sqrt_f);
};
Expand All @@ -52,24 +57,18 @@ int main() {
solver.mesh = create_mesh(mesh_name);
solver.backend = create_backend(backend_name, *solver.mesh);

// Copy initial relative pose to the body frame
// SoA_3D_t<f32> initial_pose;
// initial_pose.resize(solver.mesh->nb_vertices_wing());
// initial_pose.x = solver.mesh->v.x;
// initial_pose.y = solver.mesh->v.y;
// initial_pose.z = solver.mesh->v.z;

std::vector<f32> vec_dt; // pre-calculated timesteps
std::vector<f32> vec_t; // pre-calculated timesteps

// Pre-calculate timesteps to determine wake size
vec_t.push_back(0.0f);
for (f32 t = 0.0f; t < t_final;) {
const f32 dt = (0.5f * (solver.mesh->chord_length(0) + solver.mesh->chord_length(1))) / absolute_velocity(t, {0.f, 0.f, 0.f, 1.f});
vec_dt.push_back(dt);
t += dt;
vec_t.push_back(t);
}

// TODO: think about splitting mesh body and wake mesh
// solver.mesh->alloc(vec_dt.size()+1); // +1 for the initial pose
solver.mesh->resize_wake(vec_t.size()-1); // +1 for the initial pose

// f32 t = 0.0f;
// for (u64 i = 0; i < vec_dt.size(); i++) {
Expand Down
9 changes: 6 additions & 3 deletions vlm/include/vlm_mesh.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,23 +71,26 @@ class Mesh {
u64 nc = 1; // number of panels chordwise
u64 ns = 1; // number of panels spanwise
u64 nw; // number of wake panels chordwise (max capacity)
u64 current_nw = 0; // current number of built wake panels
u64 current_nw = 0; // current number of built wake panels (temporary)

// Analytical quanities when provided
// ---------------------
f32 s_ref = 0.0f; // reference area (of wing)
f32 c_ref = 0.0f; // reference chord (of wing)
linalg::alias::float3 ref_pt = {0.25f, 0.0f, 0.0f};

linalg::alias::float4x4 frame = linalg::identity; // transformation matrix
linalg::alias::float3 ref_pt = {0.25f, 0.0f, 0.0f}; // TODO: deprecate

void update_wake(const linalg::alias::float3& u_inf); // update wake vertices
void correction_high_aoa(f32 alpha_rad); // correct collocation points for high aoa
void compute_connectivity(); // fill offsets, connectivity
void compute_metrics_wing(); // fill colloc, normal, area
void compute_metrics_wake();
void compute_metrics_i(u64 i);
void transform(const linalg::alias::float4x4& transform);
void shed_wake();
void move(const linalg::alias::float4x4& transform);
void resize_wake(const u64 nw);

u64 nb_panels_wing() const;
u64 nb_panels_total() const;
u64 nb_vertices_wing() const;
Expand Down
19 changes: 13 additions & 6 deletions vlm/src/vlm_mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -368,25 +368,32 @@ void Mesh::io_read(const std::string& filename) {
throw std::runtime_error("Failed to open mesh file");
}
}

void Mesh::transform(const linalg::alias::float4x4& transform) {
// like a functional map but over multidimensional elements

void Mesh::move(const linalg::alias::float4x4& transform) {
for (u64 i = 0; i < nb_vertices_wing(); i++) {
linalg::alias::float4 global_vertex{v.x[i], v.y[i], v.z[i], 1.f};
linalg::alias::float4 local_vertex = frame.col(3);

const linalg::alias::float4 transformed_pt = linalg::mul(transform, linalg::alias::float4{v.x[i], v.y[i], v.z[i], 1.f});
v.x[i] = transformed_pt.x;
v.y[i] = transformed_pt.y;
v.z[i] = transformed_pt.z;
}
}

void Mesh::resize_wake(const u64 nw_) {
nw = nw_;
alloc();
}

// Duplicate the trailing edge vertices and release them in the buffer
void Mesh::shed_wake() {
const u64 src_begin = nb_vertices_wing() - ns - 1;
assert(current_nw < nw); // check if we have capacity
const u64 src_begin = nb_vertices_wing() - (ns + 1);
const u64 src_end = nb_vertices_wing();
const u64 dst_begin = (nc + 1 + current_nw) * (ns + 1);
const u64 dst_begin = (nc + nw - current_nw) * (ns + 1);
std::copy(v.x.data() + src_begin, v.x.data() + src_end, v.x.data() + dst_begin);
std::copy(v.y.data() + src_begin, v.y.data() + src_end, v.y.data() + dst_begin);
std::copy(v.z.data() + src_begin, v.z.data() + src_end, v.z.data() + dst_begin);
current_nw++;
}

2 changes: 1 addition & 1 deletion xmake.lua
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ add_rules("mode.debug", "mode.release", "mode.releasedbg", "mode.asan")
-- set_toolchains("cuda")

-- set_toolset("cxx", "clang")
-- set_policy("build.sanitizer.address", true)
set_policy("build.sanitizer.address", true)
set_policy("build.warning", true)
set_policy("build.cuda.devlink", true) -- magic
set_policy("run.autobuild", true)
Expand Down

0 comments on commit a9d8923

Please sign in to comment.