diff --git a/data/displacement_real.py b/data/displacement_real.py new file mode 100644 index 0000000..c5bf895 --- /dev/null +++ b/data/displacement_real.py @@ -0,0 +1,78 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.animation as animation +from mpl_toolkits.mplot3d.art3d import Poly3DCollection + +def coords_to_verts(coords, nc, ns): + assert coords.shape == (3, (ns+1)*(nc+1)) + verts = [] + if (nc > 0): + for ia in range(nc*ns): + v0 = ia + ia // ns + v1 = v0 + 1 + v3 = v0 + ns+1 + v2 = v3 + 1 + idx = [v0, v1, v2, v3] + verts.append([[coords[0, i], coords[1, i], coords[2, i]] for i in idx]) + + return verts + +coords_wing = [] +coords_wake = [] + +with open("build/windows/x64/release/wing_data.txt") as f: + # first line + nc, ns = map(int, f.readline().split()) + timesteps = int(f.readline()) + for i in range(timesteps): + coords = np.zeros((3, (ns+1)*(nc+1)), dtype=np.float32) + f.readline() # skip newline + coords[0, :] = np.array(list(map(float, f.readline().split()))) + coords[1, :] = np.array(list(map(float, f.readline().split()))) + coords[2, :] = np.array(list(map(float, f.readline().split()))) + coords_wing.append(coords) + +with open("build/windows/x64/release/wake_data.txt") as f: + # first line + f.readline() # skip first line + f.readline() # skip second line + for i in range(timesteps): + coords = np.zeros((3, (ns+1)*(i+1)), dtype=np.float32) + f.readline() # skip newline + coords[0, :] = np.array(list(map(float, f.readline().split()))) + coords[1, :] = np.array(list(map(float, f.readline().split()))) + coords[2, :] = np.array(list(map(float, f.readline().split()))) + coords_wake.append(coords) + +# # Setup animation +fig = plt.figure() +ax = fig.add_subplot(111, projection='3d') + +# Setting the initial view to be isometric +ax.view_init(elev=30, azim=135) # Isometric view angle +ax.set_xlabel('X axis') +ax.set_ylabel('Y axis') +ax.set_zlabel('Z axis') +ax.set_xlim(-30, 0) +ax.set_ylim(-15, 15) +ax.set_zlim(-15, 15) +ax.invert_xaxis() # Invert x axis +ax.invert_yaxis() # Invert y axis + +poly_wing = Poly3DCollection(coords_to_verts(coords_wing[0], nc, ns), facecolors=['b'], edgecolor='k') +poly_wake = Poly3DCollection(coords_to_verts(coords_wake[0], 0, ns), facecolors=['r'], edgecolor='k') +ax.add_collection3d(poly_wing) +ax.add_collection3d(poly_wake) + +def update(frame, poly_wing, poly_wake): + global coords_wing, coords_wake + + poly_wing.set_verts(coords_to_verts(coords_wing[frame], nc, ns)) + poly_wake.set_verts(coords_to_verts(coords_wake[frame], frame, ns)) + + return poly_wing, poly_wake + +ani = animation.FuncAnimation(fig, update, frames=np.arange(0, timesteps), fargs=(poly_wing, poly_wake), blit=False, repeat=False) +# ani.save('animation.mp4', fps=30, extra_args=['-vcodec', 'libx264']) + +plt.show() diff --git a/tests/test_uvlm_theodorsen.cpp b/tests/test_uvlm_theodorsen.cpp index 15f5003..effddc5 100644 --- a/tests/test_uvlm_theodorsen.cpp +++ b/tests/test_uvlm_theodorsen.cpp @@ -11,6 +11,8 @@ #include "vlm_data.hpp" #include "vlm_utils.hpp" +#define DEBUG_DISPLACEMENT_DATA + using namespace vlm; using namespace linalg::ostream_overloads; @@ -47,14 +49,24 @@ class Kinematics { } }; +template +void dump_buffer(std::ofstream& stream, T* start, T* end) { + for (T* it = start; it != end; it++) { + stream << *it << " "; + } + stream << "\n"; +} + int main() { - const std::vector meshes = {"../../../../mesh/infinite_rectangular_2x8.x"}; + // const std::vector meshes = {"../../../../mesh/infinite_rectangular_2x8.x"}; + const std::vector meshes = {"../../../../mesh/rectangular_4x4.x"}; + const std::vector backends = get_available_backends(); auto solvers = tiny::make_combination(meshes, backends); // Define simulation length - const f32 t_final = 10.0f; + const f32 t_final = 30.0f; Kinematics kinematics{}; @@ -63,7 +75,7 @@ int main() { return linalg::translation_matrix(linalg::alias::float3{ 0.0f, 0.0f, - std::sin(0.2f * t) + 2.f * std::sin(0.2f * t) }); }); @@ -76,37 +88,73 @@ int main() { }); for (const auto& [mesh_name, backend_name] : solvers) { - std::unique_ptr mesh = create_mesh(mesh_name); + const std::unique_ptr mesh = create_mesh(mesh_name); // Pre-calculate timesteps to determine wake size std::vector vec_t; // timesteps vec_t.push_back(0.0f); for (f32 t = 0.0f; t < t_final;) { - const f32 dt = (0.5f * (mesh->chord_length(0) + mesh->chord_length(1))) / kinematics.velocity_magnitude(t, {0.f, 0.f, 0.f, 1.f}); + // TODO: this is currently not accurate for rotational motion + const f32 dt = mesh->panel_length(mesh->nc-1, 0) / kinematics.velocity_magnitude(t, {0.f, 0.f, 0.f, 1.f}); t += dt; vec_t.push_back(t); } + #ifdef DEBUG_DISPLACEMENT_DATA + std::ofstream wing_data("wing_data.txt"); + std::ofstream wake_data("wake_data.txt"); + + wing_data << mesh->nc << " " << mesh->ns << "\n"; + wing_data << vec_t.size() - 1 << "\n\n"; + + wake_data << mesh->ns << "\n"; + wake_data << vec_t.size() - 1 << "\n\n"; + #endif + mesh->resize_wake(vec_t.size()-1); // +1 for the initial pose - std::unique_ptr backend = create_backend(backend_name, *mesh); // create after mesh has been resized + const std::unique_ptr backend = create_backend(backend_name, *mesh); // create after mesh has been resized // Precompute the LHS since wing geometry is constant - FlowData flow_dummy{0.0f, 0.0f, 1.0f, 1.0f}; + const FlowData flow_dummy{0.0f, 0.0f, 1.0f, 1.0f}; backend->compute_lhs(flow_dummy); backend->lu_factor(); + + // Unsteady loop + std::cout << "Timesteps: " << vec_t.size() << "\n"; for (u64 i = 0; i < vec_t.size()-1; i++) { + assert(mesh->current_nw == i); // dumb assert + + #ifdef DEBUG_DISPLACEMENT_DATA + dump_buffer(wing_data, mesh->v.x.data(), mesh->v.x.data() + mesh->nb_vertices_wing()); + dump_buffer(wing_data, mesh->v.y.data(), mesh->v.y.data() + mesh->nb_vertices_wing()); + dump_buffer(wing_data, mesh->v.z.data(), mesh->v.z.data() + mesh->nb_vertices_wing()); + wing_data << "\n"; + + const u64 wake_start = (mesh->nc + mesh->nw - i) * (mesh->ns + 1); + const u64 wake_end = mesh->nb_vertices_total(); + std::cout << "Buffer size: " << mesh->v.x.size() << " | " << wake_start << " | " << wake_end << std::endl; + + dump_buffer(wake_data, mesh->v.x.data() + wake_start, mesh->v.x.data() + wake_end); + dump_buffer(wake_data, mesh->v.y.data() + wake_start, mesh->v.y.data() + wake_end); + dump_buffer(wake_data, mesh->v.z.data() + wake_start, mesh->v.z.data() + wake_end); + wake_data << "\n"; + #endif + + const f32 t = vec_t[i]; + const f32 dt = vec_t[i+1] - t; auto base_vertex = mesh->get_v0(0); auto base_velocity = kinematics.velocity(vec_t[i], {base_vertex[0], base_vertex[1], base_vertex[2], 1.0f}); - FlowData flow{linalg::alias::float3{base_velocity[0], base_velocity[1], base_velocity[2]}, 1.0f}; + const FlowData flow{linalg::alias::float3{base_velocity[0], base_velocity[1], base_velocity[2]}, 1.0f}; backend->compute_rhs(flow); backend->add_wake_influence(flow); backend->lu_solve(); backend->compute_delta_gamma(); - // compute d gamma / dt - const f32 cl_steady = backend->compute_coefficient_cl(flow); - // const f32 cl_unsteady = backend->compute_coefficient_unsteady_cl(flow); // TODO: implement - // std::printf("t: %f, CL: %f\n", vec_t[i], cl_steady + cl_unsteady); - mesh->move(kinematics.relative_displacement(vec_t[i], vec_t[i+1])); + if (i > 0) { + // TODO: this should take a vector of local velocities magnitude because it can be different for each point on the mesh + const f32 cl_unsteady = backend->compute_coefficient_unsteady_cl(flow, dt, mesh->s_ref, 0, mesh->ns); + std::printf("t: %f, CL: %f\n", t, cl_unsteady); + } + mesh->move(kinematics.relative_displacement(t, t+dt)); backend->shed_gamma(); } } diff --git a/vlm/backends/cpu/include/vlm_backend_cpu.hpp b/vlm/backends/cpu/include/vlm_backend_cpu.hpp index 75f023a..9fa83a1 100644 --- a/vlm/backends/cpu/include/vlm_backend_cpu.hpp +++ b/vlm/backends/cpu/include/vlm_backend_cpu.hpp @@ -12,6 +12,7 @@ class BackendCPU : public Backend { std::vector rhs; std::vector ipiv; std::vector gamma; // ncw * ns + std::vector gamma_prev; // nc * ns (previous timestep gamma) std::vector delta_gamma; std::vector trefftz_buffer; @@ -26,6 +27,7 @@ class BackendCPU : public Backend { void lu_factor() override; void lu_solve() override; f32 compute_coefficient_cl(const FlowData& flow, const f32 area, const u64 j, const u64 n) override; + f32 compute_coefficient_unsteady_cl(const FlowData& flow, f32 dt, const f32 area, const u64 j, const u64 n) override; linalg::alias::float3 compute_coefficient_cm(const FlowData& flow, const f32 area, const f32 chord, const u64 j, const u64 n) override; f32 compute_coefficient_cd(const FlowData& flow, const f32 area, const u64 j, const u64 n) override; void compute_delta_gamma() override; diff --git a/vlm/backends/cpu/src/vlm_backend_cpu.cpp b/vlm/backends/cpu/src/vlm_backend_cpu.cpp index 9978fa7..d3a6096 100644 --- a/vlm/backends/cpu/src/vlm_backend_cpu.cpp +++ b/vlm/backends/cpu/src/vlm_backend_cpu.cpp @@ -29,6 +29,7 @@ BackendCPU::BackendCPU(Mesh& mesh) : Backend(mesh) { rhs.resize(mesh.nb_panels_wing()); ipiv.resize(mesh.nb_panels_wing()); gamma.resize((mesh.nc + mesh.nw) * mesh.ns); // store wake gamma as well + gamma_prev.resize(mesh.nb_panels_wing()); delta_gamma.resize(mesh.nb_panels_wing()); trefftz_buffer.resize(mesh.ns); } @@ -142,6 +143,7 @@ void BackendCPU::shed_gamma() { Mesh& m = mesh; const u64 wake_row_start = (m.nc + m.nw - m.current_nw - 1) * m.ns; + std::copy(gamma.data(), gamma.data() + m.nb_panels_wing(), gamma_prev.data()); // store current timestep for delta_gamma std::copy(gamma.data() + m.ns * (m.nc-1), gamma.data() + m.nb_panels_wing(), gamma.data() + wake_row_start); } @@ -159,7 +161,7 @@ void BackendCPU::compute_rhs(const FlowData& flow, const std::vector& secti } void BackendCPU::lu_factor() { - tiny::ScopedTimer timer("Factor"); + const tiny::ScopedTimer timer("Factor"); const int32_t n = static_cast(mesh.nb_panels_wing()); LAPACKE_sgetrf(LAPACK_COL_MAJOR, n, n, lhs.data(), n, ipiv.data()); } @@ -172,8 +174,9 @@ void BackendCPU::lu_solve() { LAPACKE_sgetrs(LAPACK_COL_MAJOR, 'N', n, 1, lhs.data(), n, ipiv.data(), gamma.data(), n); } -f32 BackendCPU::compute_coefficient_cl(const FlowData& flow, const f32 area, - const u64 j, const u64 n) { +f32 BackendCPU::compute_coefficient_cl(const FlowData& flow, const f32 area, const u64 j, const u64 n) { + const tiny::ScopedTimer timer("Compute CL"); + assert(n > 0); assert(j >= 0 && j+n <= mesh.ns); @@ -184,15 +187,45 @@ f32 BackendCPU::compute_coefficient_cl(const FlowData& flow, const f32 area, const u64 li = u * mesh.ns + v; // linear index const linalg::alias::float3 v0 = mesh.get_v0(li); const linalg::alias::float3 v1 = mesh.get_v1(li); - const linalg::alias::float3 v3 = mesh.get_v3(li); + // const linalg::alias::float3 v3 = mesh.get_v3(li); // Leading edge vector pointing outward from wing root linalg::alias::float3 dl = v1 - v0; - const linalg::alias::float3 local_left_chord = linalg::normalize(v3 - v0); - const linalg::alias::float3 projected_vector = linalg::dot(dl, local_left_chord) * local_left_chord; - dl -= projected_vector; // orthogonal leading edge vector + // const linalg::alias::float3 local_left_chord = linalg::normalize(v3 - v0); + // const linalg::alias::float3 projected_vector = linalg::dot(dl, local_left_chord) * local_left_chord; + // dl -= projected_vector; // orthogonal leading edge vector // Distance from the center of leading edge to the reference point - const linalg::alias::float3 force = linalg::cross(flow.stream_axis, dl) * flow.rho * delta_gamma[li]; // * local flow magnitude - cl += linalg::dot(force, flow.lift_axis); + const linalg::alias::float3 force = linalg::cross(flow.freestream, dl) * flow.rho * delta_gamma[li]; + cl += linalg::dot(force, flow.lift_axis); // projection on the body lift axis + } + } + cl /= 0.5f * flow.rho * linalg::length2(flow.freestream) * area; + + return cl; +} + +f32 BackendCPU::compute_coefficient_unsteady_cl(const FlowData& flow, f32 dt, const f32 area, const u64 j, const u64 n) { + assert(n > 0); + assert(j >= 0 && j+n <= mesh.ns); + + f32 cl = 0.0f; + + for (u64 u = 0; u < mesh.nc; u++) { + for (u64 v = j; v < j + n; v++) { + const u64 li = u * mesh.ns + v; // linear index + + // Steady part: + const linalg::alias::float3 v0 = mesh.get_v0(li); + const linalg::alias::float3 v1 = mesh.get_v1(li); + // Leading edge vector pointing outward from wing root + const linalg::alias::float3 dl = v1 - v0; + const linalg::alias::float3 force_steady = flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, dl); + cl += linalg::dot(force_steady, flow.lift_axis); + + // Unsteady part (Simpson method) + const f32 gamma_dt = (gamma[li] - gamma_prev[li]) / dt; // backward difference + const linalg::alias::float3 normal{mesh.normal.x[li], mesh.normal.y[li], mesh.normal.z[li]}; + const linalg::alias::float3 force_unsteady = flow.rho * gamma_dt * mesh.area[li] * normal; + cl += linalg::dot(force_unsteady, flow.lift_axis); } } cl /= 0.5f * flow.rho * linalg::length2(flow.freestream) * area; diff --git a/vlm/include/vlm_backend.hpp b/vlm/include/vlm_backend.hpp index 31f4cf5..3ceb0a2 100644 --- a/vlm/include/vlm_backend.hpp +++ b/vlm/include/vlm_backend.hpp @@ -19,6 +19,7 @@ class Backend { virtual void lu_factor() = 0; virtual void lu_solve() = 0; virtual f32 compute_coefficient_cl(const FlowData& flow, const f32 area, const u64 j, const u64 n) = 0; + virtual f32 compute_coefficient_unsteady_cl(const FlowData& flow, f32 dt, const f32 area, const u64 j, const u64 n) = 0; f32 compute_coefficient_cl(const FlowData& flow); virtual linalg::alias::float3 compute_coefficient_cm(const FlowData& flow, const f32 area, const f32 chord, const u64 j, const u64 n) = 0; linalg::alias::float3 compute_coefficient_cm(const FlowData& flow); diff --git a/vlm/include/vlm_mesh.hpp b/vlm/include/vlm_mesh.hpp index 1cc7f3a..3b715b8 100644 --- a/vlm/include/vlm_mesh.hpp +++ b/vlm/include/vlm_mesh.hpp @@ -97,6 +97,7 @@ class Mesh { f32 panels_area(const u64 i, const u64 j, const u64 m, const u64 n) const; f32 panels_area_xy(const u64 i, const u64 j, const u64 m, const u64 n) const; f32 panel_width_y(const u64 i, const u64 j) const; + f32 panel_length(const u64 i, const u64 j) const; f32 strip_width(const u64 j) const; f32 chord_length(const u64 j) const; // vertex idx f32 chord_mean(const u64 j, const u64 n) const; diff --git a/vlm/src/vlm_mesh.cpp b/vlm/src/vlm_mesh.cpp index 2f3b4f6..47a2358 100644 --- a/vlm/src/vlm_mesh.cpp +++ b/vlm/src/vlm_mesh.cpp @@ -138,14 +138,25 @@ f32 Mesh::panels_area_xy(const u64 i, const u64 j, const u64 m, const u64 n) con return total_area; } +f32 Mesh::panel_length(const u64 i, const u64 j) const { + // Since chordwise segments are always parallel, we can simply measure the width of the first panel + assert(i < nc); + assert(j < ns); + const linalg::alias::float3 v0 = get_v0(j + i * ns); + const linalg::alias::float3 v1 = get_v1(j + i * ns); + const linalg::alias::float3 v2 = get_v2(j + i * ns); + const linalg::alias::float3 v3 = get_v3(j + i * ns); + return 0.5f * (linalg::length(v3-v0) + linalg::length(v2-v1)); +} + /// @brief Computes the width of a single panel in pure y direction /// @param i panel index chordwise /// @param j panel index spanwise f32 Mesh::panel_width_y(const u64 i, const u64 j) const { - // Since chordwise segments are always parallel, we can simply measure the width of the first panel assert(i < nc); assert(j < ns); const u64 ld = ns + 1; + return v.y[j + 1 + i * ld] - v.y[j + i * ld]; }