more debugging

config/swept.vlm

@@ -5,9 +5,10 @@
 # Y = 10
 
 <solver>
-ref_point = [0.25, 0.0, 0.0]
 alphas = [5]
 backend = cpu
+sigma_vatistas = 0.0
 
 <mesh>
+ref_pt = [0.25, 0.0, 0.0]
 wake_included = false

data/theodorsen.py

@@ -39,8 +39,9 @@ def atime(t: float): return 2. * u_inf * t / c
 amplitudes = [0.1, 0.1, 0.1] 
 reduced_frequencies = [0.5, 0.75, 1.5]
 
-t_final = 40
-t = np.linspace(0, t_final, 500)
+t_final = 30
+nb_pts = 500
+t = np.linspace(0, t_final, nb_pts)
 
 fig, axs = plt.subplot_mosaic(
     [["time"],["heave"]],  # Disposition des graphiques
@@ -87,21 +88,23 @@ def cl_theodorsen(t: float): # using Wagner functions and Kholodar formulation
     coord_z = np.array([heave(ti) / (2*b) for ti in t])
 
     axs["time"].plot(t, cl, label=f"k={k}")
-    axs["heave"].plot(coord_z[len(cl)//2:], cl[len(cl)//2:], label=f"k={k}")
+    axs["heave"].plot(coord_z[int(nb_pts//2):], cl[int(nb_pts//2):], label=f"k={k}")
 
 uvlm_cl = []
 uvlm_t = []
 uvlm_z = []
-with open("build/windows/x64/release/cl_data.txt", "r") as f:
+with open("build/windows/x64/debug/cl_data.txt", "r") as f:
     for line in f:
-        t, z, cl = line.split()
-        uvlm_t.append(float(t))
+        time_step, z, cl = line.split()
+        uvlm_t.append(float(time_step))
         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["heave"].plot(uvlm_z[len(uvlm_cl)//2:], uvlm_cl[len(uvlm_cl)//2:], label="UVLM (k=0.5)", linestyle="--")
+axs["time"].plot(uvlm_t, uvlm_z, label="Heave (k=0.5)", linestyle="--")
+axs["heave"].plot(uvlm_z[len(uvlm_cl)//4:], uvlm_cl[len(uvlm_cl)//4:], label="UVLM (k=0.5)", linestyle="--")
+
+# axs["time"].plot(t, [0.548311] * len(t), label="VLM (alpha=5)")
 
 axs["time"].set_xlabel('t')
 axs["time"].set_ylabel('CL')

headeronly/linalg.h

@@ -531,9 +531,21 @@ namespace linalg
     enum z_range { neg_one_to_one, zero_to_one };   // Should projection matrices map z into the range of [-1,1] or [0,1]?
     template<class T> vec<T,4>   rotation_quat     (const vec<T,3> & axis, T angle)         { return {axis*std::sin(angle/2), std::cos(angle/2)}; }
     template<class T> vec<T,4>   rotation_quat     (const mat<T,3,3> & m);
+    template<class T> mat<T,3,3> skew_matrix       (const vec<T,3> & a)                      { return {{0, a.z, -a.y}, {-a.z, 0, a.x}, {a.y, -a.x, 0}}; }
     template<class T> mat<T,4,4> translation_matrix(const vec<T,3> & translation)           { return {{1,0,0,0},{0,1,0,0},{0,0,1,0},{translation,1}}; }
     template<class T> mat<T,4,4> rotation_matrix   (const vec<T,4> & rotation)              { return {{qxdir(rotation),0}, {qydir(rotation),0}, {qzdir(rotation),0}, {0,0,0,1}}; }
-    template<class T> mat<T,4,4> rotation_matrix   (const vec<T,3> & axis, const vec<T,3> & p, T angle) {}
+    template<class T> mat<T,4,4> rotation_matrix   (const vec<T,3> & point, const vec<T,3> & axis, T angle) {
+        //mat<T,4,4> matrix = identity;
+        const mat<T,3,3> skew_mat = skew_matrix(axis);
+        const mat<T,3,3> rodrigues = identity + std::sin(angle)*skew_mat + (1-std::cos(angle))*mul(skew_mat, skew_mat);
+        const vec<T,3> trans_part = mul((identity - rodrigues), point);
+        return {
+            {rodrigues.x.x, rodrigues.x.y, rodrigues.x.z, 0}, // 1st col
+            {rodrigues.y.x, rodrigues.y.y, rodrigues.y.z, 0},
+            {rodrigues.z.x, rodrigues.z.y, rodrigues.z.z, 0},
+            {trans_part.x, trans_part.y, trans_part.z, 1}
+        };
+    }
     template<class T> mat<T,4,4> scaling_matrix    (const vec<T,3> & scaling)               { return {{scaling.x,0,0,0}, {0,scaling.y,0,0}, {0,0,scaling.z,0}, {0,0,0,1}}; }
     template<class T> mat<T,4,4> pose_matrix       (const vec<T,4> & q, const vec<T,3> & p) { return {{qxdir(q),0}, {qydir(q),0}, {qzdir(q),0}, {p,1}}; }
     template<class T> mat<T,4,4> lookat_matrix     (const vec<T,3> & eye, const vec<T,3> & center, const vec<T,3> & view_y_dir, fwd_axis fwd = neg_z);

tests/test_uvlm_theodorsen.cpp

@@ -83,27 +83,30 @@ int main() {
 
     Kinematics kinematics{};
 
-    // Define wing motion
+    // Periodic heaving
     kinematics.add([=](f32 t) {
         return linalg::translation_matrix(linalg::alias::float3{
+            -u_inf*t, // freestream
             0.0f,
-            0.0f,
-            amplitude * std::sin(omega* t)
+            amplitude * std::sin(omega* t) // heaving
         });
     });
 
-    kinematics.add([=](f32 t) {
-        return linalg::translation_matrix(linalg::alias::float3{
-            -u_inf*t,
-            0.0f,
-            0.0f
-        });
-    });
+    // Sudden acceleration
+    // const f32 alpha = to_radians(5.0f);
+    // kinematics.add([=](f32 t) {
+    //     return linalg::translation_matrix(linalg::alias::float3{
+    //         -u_inf*cos(alpha)*t,
+    //         0.0f,
+    //         -u_inf*sin(alpha)*t
+    //     });
+    // });
 
     for (const auto& [mesh_name, backend_name] : solvers) {
         const std::unique_ptr<Mesh> mesh = create_mesh(mesh_name);
 
         // Pre-calculate timesteps to determine wake size
+        // Note: calculation should be made on the four corners of the wing
         std::vector<f32> vec_t; // timesteps
         vec_t.push_back(0.0f);
         for (f32 t = 0.0f; t < t_final;) {
@@ -136,8 +139,6 @@ int main() {
         // 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());
@@ -179,10 +180,10 @@ int main() {
                 cl_data << t << " " << mesh->v.z[0] << " " << cl_unsteady << "\n";
                 #endif
             }
+            backend->shed_gamma(); // shed before moving & incrementing currentnw
             mesh->move(kinematics.relative_displacement(t, t+dt));
-            backend->shed_gamma();
         }
     }
-
+  
     return 0;
 }

vlm/backends/cpu/src/vlm_backend_cpu.cpp

@@ -19,6 +19,7 @@
 #include <cblas.h>
 
 using namespace vlm;
+using namespace linalg::ostream_overloads;
 
 BackendCPU::~BackendCPU() = default; // Destructor definition
 
@@ -125,49 +126,49 @@ void BackendCPU::add_wake_influence(const FlowData& flow) {
     };
 
     // loop over wake rows
-    // for (u64 i = 0; i < mesh.current_nw; i++) {
-    //     const u64 wake_row_start = (m.nc + m.nw - i - 1) * m.ns;
-    //     std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out 
-    //     // Actually fill the wake buffer
-    //     // parallel for
-    //     for (u64 j = 0; j < mesh.ns; j++) { // loop over columns
-    //         const u64 lidx = wake_row_start + j; // TODO: replace this ASAP
-    //         ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas);
-    //     }
+    for (u64 i = 0; i < mesh.current_nw; i++) {
+        const u64 wake_row_start = (m.nc + m.nw - i - 1) * m.ns;
+        std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out 
+        // Actually fill the wake buffer
+        // parallel for
+        for (u64 j = 0; j < mesh.ns; j++) { // loop over columns
+            const u64 lidx = wake_row_start + j; // TODO: replace this ASAP
+            ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas);
+        }
 
-    //     cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, -1.0f, wake_buffer.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, rhs.data(), 1);
-    // }
+        cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, -1.0f, wake_buffer.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, rhs.data(), 1);
+    }
 
-    u64 idx = 0;
-    u64 wake_row_start = (m.nc + m.nw - 1) * m.ns;
+    // u64 idx = 0;
+    // u64 wake_row_start = (m.nc + m.nw - 1) * m.ns;
 
-    auto init = taskflow.placeholder();
-    auto cond = taskflow.emplace([&]{
-        return idx < mesh.current_nw ? 0 : 1; // 0 means continue, 1 means break
-    });
-    auto zero_buffer = taskflow.emplace([&]{
-        std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out 
-    });
-    auto wake_influence = taskflow.for_each_index((u64)0, m.ns, [&] (u64 j) {
-        const u64 lidx = wake_row_start + j;
-        ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas);
-    });
-    auto back = taskflow.emplace([&]{
-        cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, -1.0f, wake_buffer.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, rhs.data(), 1);
+    // auto init = taskflow.placeholder();
+    // auto cond = taskflow.emplace([&]{
+    //     return idx < mesh.current_nw ? 0 : 1; // 0 means continue, 1 means break
+    // });
+    // auto zero_buffer = taskflow.emplace([&]{
+    //     std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out 
+    // });
+    // auto wake_influence = taskflow.for_each_index((u64)0, m.ns, [&] (u64 j) {
+    //     const u64 lidx = wake_row_start + j;
+    //     ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas);
+    // });
+    // auto back = taskflow.emplace([&]{
+    //     cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, -1.0f, wake_buffer.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, rhs.data(), 1);
 
-        idx++;
-        wake_row_start -= m.ns;
-        return 0; // 0 means continue
-    });
-    auto sync = taskflow.placeholder();
+    //     idx++;
+    //     wake_row_start -= m.ns;
+    //     return 0; // 0 means continue
+    // });
+    // auto sync = taskflow.placeholder();
 
-    init.precede(cond);
-    cond.precede(zero_buffer, sync);
-    zero_buffer.precede(wake_influence);
-    wake_influence.precede(back);
-    back.precede(cond);
+    // init.precede(cond);
+    // cond.precede(zero_buffer, sync);
+    // zero_buffer.precede(wake_influence);
+    // wake_influence.precede(back);
+    // back.precede(cond);
 
-    Executor::get().run(taskflow).wait();
+    // Executor::get().run(taskflow).wait();
 }
 
 void BackendCPU::shed_gamma() {
@@ -250,20 +251,20 @@ f32 BackendCPU::compute_coefficient_unsteady_cl(const FlowData& flow, f32 dt, co
             const linalg::alias::float3 v2 = mesh.get_v2(li);
             const linalg::alias::float3 v3 = mesh.get_v3(li);
 
-            linalg::alias::float3 force_steady = {0.0f, 0.0f, 0.0f};
-            force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v1 - v0);
-            force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v2 - v1);
-            force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v3 - v2);
-            force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v0 - v3);
+            linalg::alias::float3 force = {0.0f, 0.0f, 0.0f};
+            force += flow.rho * gamma[li] * linalg::cross(flow.freestream, linalg::normalize(v1 - v0));
+            force += flow.rho * gamma[li] * linalg::cross(flow.freestream, linalg::normalize(v2 - v1));
+            force += flow.rho * gamma[li] * linalg::cross(flow.freestream, linalg::normalize(v3 - v2));
+            force += flow.rho * gamma[li] * linalg::cross(flow.freestream, linalg::normalize(v0 - v3));
 
             // Leading edge vector pointing outward from wing root
-            cl += linalg::dot(force_steady, flow.lift_axis);
+            //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);
+            force += flow.rho * gamma_dt * mesh.area[li] * normal;
+            cl += linalg::dot(force, flow.lift_axis);
         }
     }
     cl /= 0.5f * flow.rho * linalg::length2(flow.freestream) * area;