debug efforts

data/displacement.py

@@ -17,15 +17,14 @@ def displacement(self, t):
 
     def relative_displacement(self, t0, t1):
         return self.displacement(t1) @ np.linalg.inv(self.displacement(t0))
-
-    # f32 absolute_velocity(f32 t, const linalg::alias::float4& vertex) {
-    #     return linalg::length((linalg::mul(displacement(t+EPS_sqrt_f), vertex)-linalg::mul(displacement(t), vertex))/EPS_sqrt_f);
-    # }
 
-    def absolute_velocity(self, t, vertex):
+    def velocity(self, t, vertex):
         EPS_sqrt_f = np.sqrt(1.19209e-07)
-        return np.linalg.norm((self.relative_displacement(t, t+EPS_sqrt_f) @ vertex - vertex) / EPS_sqrt_f)
+        return(self.relative_displacement(t, t+EPS_sqrt_f) @ vertex - vertex) / EPS_sqrt_f
 
+    def absolute_velocity(self, t, vertex):
+        return np.linalg.norm(self.velocity(t, vertex))
+
 def skew_matrix(v):
     return np.array([
         [0, -v[2], v[1]],
@@ -48,21 +47,11 @@ def rotation_matrix(center, axis, theta):
     mat[:3, 3] = (np.identity(3) - rodrigues) @ center # matvec
     return mat
 
-def rotation_matrix2(axis, theta):
-    axis = np.array(axis)
-    axis = axis / np.sqrt(np.dot(axis, axis))
-    a = np.cos(theta / 2.0)
-    b, c, d = -axis * np.sin(theta / 2.0)
-    return np.array([[a*a + b*b - c*c - d*d, 2*(b*c - a*d), 2*(b*d + a*c), 0],
-                     [2*(b*c + a*d), a*a + c*c - b*b - d*d, 2*(c*d - a*b), 0],
-                     [2*(b*d - a*c), 2*(c*d + a*b), a*a + d*d - b*b - c*c, 0],
-                     [0, 0, 0, 1]])
-
 def move_vertices(vertices, displacement):
     return displacement @ vertices
 
 def displacement_wing(t): return translation_matrix([0, 0, np.sin(0.9 * t)])
-def displacement_freestream(t): return translation_matrix([-2 * t, 0, 0])
+def displacement_freestream(t): return translation_matrix([-1 * t, 0, 0])
 # def displacement_rotor(t, frame): 
 #     return rotation_matrix(frame @ [0, 0, 0, 1], frame @ [0, 0, 1, 0], 1 * t)
 def displacement_rotor(t): 
@@ -73,9 +62,9 @@ def pitching(t):
 kinematics = Kinematics()
 kinematics.add_joint(displacement_freestream, np.identity(4))
 # kinematics.add_joint(displacement_wing, np.identity(4))
-kinematics.add_joint(pitching, np.identity(4))
+# kinematics.add_joint(pitching, np.identity(4))
 
-dt = 0.1
+dt = 0.2
 t_final = 15
 
 # vertices of a single panel (clockwise) (initial position) (global coordinates)
@@ -106,33 +95,36 @@ def pitching(t):
 line, = ax.plot3D(vertices[0, :], vertices[1, :], vertices[2, :], '-') 
 scatter = ax.scatter(vertices[0, :], vertices[1, :], vertices[2, :], c='r', marker='o')
 
-current_frame = 0
-def update(frame):
-    global vertices, kinematics, current_frame
-    t = frame * dt
+print(kinematics.velocity(0.2, [0, 0, 0, 1]))
+
+# current_frame = 0
+# def update(frame):
+#     global vertices, kinematics, current_frame
+#     t = frame * dt
 
-    vertices_velocity = np.array([kinematics.absolute_velocity(t, vertex) for vertex in vertices.T])
-    if frame == current_frame: # otherwise invalid velocity value
-        print(f"frame: {frame} | vel: {vertices_velocity[:-1]}")
+#     vertices_velocity = np.array([kinematics.absolute_velocity(t, vertex) for vertex in vertices.T])
+#     if frame == current_frame: # otherwise invalid velocity value
+#         print(f"velocity: {kinematics.velocity(t, vertices[:, 0])}")
+#         print(f"frame: {frame} | vel: {vertices_velocity[:-1]}")
 
-    norm = plt.Normalize(vertices_velocity.min(), vertices_velocity.max())
-    colors = cm.viridis(norm(vertices_velocity))
+#     norm = plt.Normalize(vertices_velocity.min(), vertices_velocity.max())
+#     colors = cm.viridis(norm(vertices_velocity))
 
-    # Update the line object for 3D
-    line.set_data(vertices[0, :], vertices[1, :])  # y and z for 2D part of set_data
-    line.set_3d_properties(vertices[2, :])  # x for the 3rd dimension
+#     # Update the line object for 3D
+#     line.set_data(vertices[0, :], vertices[1, :])  # y and z for 2D part of set_data
+#     line.set_3d_properties(vertices[2, :])  # x for the 3rd dimension
 
-    scatter._offsets3d = (vertices[0, :], vertices[1, :], vertices[2, :])
-    scatter.set_facecolor(colors)
+#     scatter._offsets3d = (vertices[0, :], vertices[1, :], vertices[2, :])
+#     scatter.set_facecolor(colors)
 
-    if (frame == current_frame): # fix double frame 0 issue
-        print(f"t = {t:.2f}/{t_final}", end='\r')
-        vertices = move_vertices(vertices, kinematics.relative_displacement(t, t+dt))
-        current_frame += 1
+#     if (frame == current_frame): # fix double frame 0 issue
+#         print(f"t = {t:.2f}/{t_final}", end='\r')
+#         vertices = move_vertices(vertices, kinematics.relative_displacement(t, t+dt))
+#         current_frame += 1
 
-    return line, scatter
+#     return line, scatter
 
-ani = animation.FuncAnimation(fig, update, frames=np.arange(0, t_final/dt), blit=False, repeat=False)
-# ani.save('animation.mp4', fps=30, extra_args=['-vcodec', 'libx264'])
+# ani = animation.FuncAnimation(fig, update, frames=np.arange(0, t_final/dt), blit=False, repeat=False)
+# # ani.save('animation.mp4', fps=30, extra_args=['-vcodec', 'libx264'])
 
-plt.show()
+# plt.show()

data/theodorsen.py

@@ -14,6 +14,7 @@ def derivative2(f, x):
 def solve_ivp(x0: float, s0: float, sf: float, f: callable):
     return spi.solve_ivp(f, [s0, sf], [x0]).y[-1] # return only the result at t=sf
 
+# Some info in Katz Plotkin p414 (eq 13.73a)
 # Jone's approximation of Wagner function
 b0 = 1
 b1 = -0.165
@@ -36,7 +37,7 @@ def solve_ivp(x0: float, s0: float, sf: float, f: callable):
 fig, axs = plt.subplot_mosaic(
     [["time"],["heave"]],  # Disposition des graphiques
     constrained_layout=True,  # Demander à Matplotlib d'essayer d'optimiser la disposition des graphiques pour que les axes ne se superposent pas
-    figsize=(16, 9),  # Ajuster la taille de la figure (x,y)
+    figsize=(11, 6),  # Ajuster la taille de la figure (x,y)
 )
 
 for amp, k in zip(amplitudes, reduced_frequencies):
@@ -66,14 +67,18 @@ def cl_theodorsen(t: float):
 
 uvlm_cl = []
 uvlm_t = []
+uvlm_z = []
 with open("build/windows/x64/release/cl_data.txt", "r") as f:
     for line in f:
-        t, cl = line.split()
+        t, z, cl = line.split()
         uvlm_t.append(float(t))
+        uvlm_z.append(float(z))
         uvlm_cl.append(float(cl))
 
 axs["time"].plot(uvlm_t, uvlm_cl, label="UVLM (k=0.5)", linestyle="--")
 
+axs["time"].plot(uvlm_t, uvlm_z, label="UVLM z (k=0.5)", linestyle="--")
+
 axs["time"].set_xlabel('t')
 axs["time"].set_ylabel('CL')
 axs["time"].grid(True, which='both', linestyle=':', linewidth=1.0, color='gray')

headeronly/linalg.h

@@ -275,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 &         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; }
+        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; }
 
         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); }
@@ -293,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 &         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; }
+        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; }
 
         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); }

tests/test_uvlm_theodorsen.cpp

@@ -9,14 +9,15 @@
 
 #include "vlm.hpp"
 #include "vlm_data.hpp"
+#include "vlm_types.hpp"
 #include "vlm_utils.hpp"
 
 #define DEBUG_DISPLACEMENT_DATA
 
 using namespace vlm;
 using namespace linalg::ostream_overloads;
 
-using tmatrix = linalg::alias::float4x4; // transformation matrix
+using tmatrix = linalg::alias::float4x4; // transformation matri
 
 class Kinematics {
     public:
@@ -37,11 +38,15 @@ class Kinematics {
     }
 
     tmatrix relative_displacement(f32 t0, f32 t1) {
+        //std::printf("t0: %.10f, t1: %.10f ", t0, t1);
         return linalg::mul(displacement(t1), linalg::inverse(displacement(t0)));
     }
 
     linalg::alias::float4 velocity(f32 t, const linalg::alias::float4& vertex) {
-        return (linalg::mul(relative_displacement(t, t+EPS_sqrt_f), vertex)-vertex)/EPS_sqrt_f;
+        //const f32 EPS = std::max(10.f * t *EPS_f, EPS_f); // adaptive epsilon
+        const f32 EPS = std::max(std::sqrt(t) * EPS_sqrt_f, EPS_f);
+        return (linalg::mul(relative_displacement(t, t+EPS), vertex)-vertex)/EPS;
+        //return (linalg::mul(relative_displacement(t, t+EPS), vertex) - linalg::mul(relative_displacement(t, t-EPS), vertex))/ (2*EPS); // central diff
     }
 
     f32 velocity_magnitude(f32 t, const linalg::alias::float4& vertex) {
@@ -142,7 +147,7 @@ int main() {
 
             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;
+            // 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);
@@ -152,19 +157,25 @@ int main() {
 
             const f32 t = vec_t[i];
             const f32 dt = vec_t[i+1] - t;
+            std::cout << "\n----------------\n" << "T = " << t << "\n";
+
             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});
-            const FlowData flow{linalg::alias::float3{base_velocity[0], base_velocity[1], base_velocity[2]}, 1.0f};
+            auto base_velocity = kinematics.velocity(t, {base_vertex[0], base_vertex[1], base_vertex[2], 1.0f});
+            std::cout << "velocity: " << base_velocity << "\n";
+            std::cout << "position: " << base_vertex << "\n";
+
+            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();
             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);
+                // const f32 cl_unsteady = backend->compute_coefficient_cl(flow);
                 std::printf("t: %f, CL: %f\n", t, cl_unsteady);
                 #ifdef DEBUG_DISPLACEMENT_DATA
-                cl_data << t << " " << cl_unsteady << "\n";
+                cl_data << t << " " << mesh->v.z[0] << " " << cl_unsteady << "\n";
                 #endif
             }
             mesh->move(kinematics.relative_displacement(t, t+dt));

vlm/backends/cpu/src/vlm_backend_cpu.cpp

@@ -104,15 +104,15 @@ void kernel_cpu_rhs(u64 n, const float normal_x[], const float normal_y[], const
 }
 
 void BackendCPU::compute_rhs(const FlowData& flow) {
-    tiny::ScopedTimer timer("RHS");
+    const tiny::ScopedTimer timer("RHS");
     const Mesh& m = mesh;
 
     kernel_cpu_rhs(m.nb_panels_wing(), m.normal.x.data(), m.normal.y.data(), m.normal.z.data(), flow.freestream.x, flow.freestream.y, flow.freestream.z, rhs.data());
 }
 
 // TODO: consider changing FlowData to SolverData
 void BackendCPU::add_wake_influence(const FlowData& flow) {
-    tiny::ScopedTimer timer("Wake Influence");
+    const tiny::ScopedTimer timer("Wake Influence");
 
     tf::Taskflow taskflow;
 
@@ -189,7 +189,7 @@ f32 BackendCPU::compute_coefficient_cl(const FlowData& flow, const f32 area, con
             const linalg::alias::float3 v1 = mesh.get_v1(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 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

xmake.lua

@@ -15,7 +15,7 @@ set_policy("run.autobuild", true)
 
 set_warnings("all")
 set_languages("c++17", "c99")
-set_runtimes("MD") -- msvc runtime library (MD/MT/MDd/MTd)
+set_runtimes("MT") -- msvc runtime library (MD/MT/MDd/MTd)
 
 -- Define backends and helper functions
 backends = {"cuda", "cpu"}