solid progress

- visualization of the displacement with displacement_real.py
- Added computation of the forces using Simpson method

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()

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<typename T>
+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<std::string> meshes = {"../../../../mesh/infinite_rectangular_2x8.x"};
+    // const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_2x8.x"};
+    const std::vector<std::string> meshes = {"../../../../mesh/rectangular_4x4.x"};
+
     const std::vector<std::string> 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> mesh = create_mesh(mesh_name);
+        const std::unique_ptr<Mesh> mesh = create_mesh(mesh_name);
 
         // Pre-calculate timesteps to determine wake size
         std::vector<f32> 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> backend = create_backend(backend_name, *mesh); // create after mesh has been resized
+        const std::unique_ptr<Backend> 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();
         }
     }

vlm/backends/cpu/include/vlm_backend_cpu.hpp

@@ -12,6 +12,7 @@ class BackendCPU : public Backend {
         std::vector<f32> rhs;
         std::vector<i32> ipiv;
         std::vector<f32> gamma; // ncw * ns
+        std::vector<f32> gamma_prev; // nc * ns (previous timestep gamma)
         std::vector<f32> delta_gamma;
         std::vector<f32> 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;

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<f32>& secti
 }
 
 void BackendCPU::lu_factor() {
-    tiny::ScopedTimer timer("Factor");
+    const tiny::ScopedTimer timer("Factor");
     const int32_t n = static_cast<int32_t>(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;

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);

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;

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];
 }