UVLM 3DOF works

data/3dof.py

@@ -186,7 +186,7 @@ def yn_func(self, y_: np.ndarray):
         yn = np.zeros_like(y)
         gamma = 0 # Cubic factor (0 for linear spring)
         yn[0, :] = self.v.mu * y[3, :] * self.v.sigma**2 * (1 + gamma * y[3, :]**2)
-        yn[2, :] = self.v.mu * ((self.v.omega_beta / self.v.omega_alpha)**2 * self.v.r_beta**2) * alpha_linear(y[5, :])      
+        yn[2, :] = self.v.mu * ((self.v.omega_beta / self.v.omega_alpha)**2 * self.v.r_beta**2) * alpha_freeplay(y[5, :])      
         return yn.reshape(y_.shape)
 
     def uncoupled_system(self, t, y: np.ndarray):
@@ -368,7 +368,7 @@ def format_subplot(fig, row, col, xlabel, ylabel):
 
     beta_peaks = []
     beta_vel = []
-    coupled_sim = False
+    coupled_sim = True
 
     for U in tqdm(U_vec): 
         # Conner
@@ -402,7 +402,7 @@ def format_subplot(fig, row, col, xlabel, ylabel):
         v.mu = v.m / (np.pi * v.rho * v.b**2)
 
         dt = 0.09
-        t_final = 3 * v.omega_alpha
+        t_final = 10 * v.omega_alpha
         vec_t = np.arange(0, t_final, dt)
         n = len(vec_t)
 
@@ -459,7 +459,7 @@ def format_subplot(fig, row, col, xlabel, ylabel):
             add_data_and_psd(fig, vec_t, np.degrees(mono.y[2, :]), "Theodorsen", 5, 2) # dbeta
             add_data_and_psd(fig, vec_t, aero_forces[2, :]/v.mu, "Theodorsen", 5, 3) # force beta
 
-            uvlm_file = Path("build/windows/x64/debug/3dof.txt")
+            uvlm_file = Path("build/windows/x64/release/3dof.txt")
             if uvlm_file.exists():
                 with open(uvlm_file, "r") as f:
                     uvlm_tsteps = int(f.readline())

tests/uvlm_2dof_modified.cpp

@@ -141,7 +141,6 @@ class UVLM_2DOF final: public Simulation {
 };
 
 UVLM_2DOF::UVLM_2DOF(const std::string& backend_name,const Assembly& assembly) : Simulation(backend_name, assembly.mesh_filenames()), m_assembly(assembly) {
-    assert(meshes.size() == 1); // single body only
     solver = backend->create_lu_solver();
     accel_solver = backend->create_lu_solver();
     alloc_buffers();
@@ -306,6 +305,13 @@ void UVLM_2DOF::update_wake_and_gamma(i64 iteration) {
         gamma_wing_i.slice(All, -1).to(gamma_wake_i.slice(All, -1-iteration));
 
         // Compute delta
+        if (m == 1) { // gross hack for missing connectivity on the flap trailing edge
+            backend->blas->axpy(
+                -1.0f,
+                gamma_wing.views()[0].slice(All, -1),
+                gamma_wing_delta_i.slice(All, 0)
+            );
+        }
         backend->blas->axpy(
             -1.0f,
             gamma_wing_i.slice(All, Range{0, -2}),
@@ -712,14 +718,24 @@ void UVLM_2DOF::run(const UVLM_2DOF_Vars& vars, f32 t_final_nd) {
 
             rhs.view().fill(0.f);
             backend->rhs_assemble_velocities(rhs.view(), normals_d.views(), velocities.views());
-            // backend->rhs_assemble_wake_influence(rhs.view(), gamma_wake.views(), colloc_d.views(), normals_d.views(), verts_wake.views(), m_assembly.lifting(), i+1);
+            // DEBUG
+            if (i == 0 || i == 1) {
+                std::cout << "After velocities\n";
+                std::cout << rhs.view() << "\n";
+
+                // for (i64 i = 0; i < m_assembly.num_surfaces(); i++) {
+                //     std::cout << velocities.views()[i] << "\n";
+                // }
+            }
+            backend->rhs_assemble_wake_influence(rhs.view(), gamma_wake.views(), colloc_d.views(), normals_d.views(), verts_wake.views(), m_assembly.lifting(), i+1);
 
             // DEBUG
             if (i == 0 || i == 1) {
-                for (i64 i = 0; i < m_assembly.num_surfaces(); i++) {
-                    // std::cout << gamma_wing.views()[i] << "\n";
-                    std::cout << velocities.views()[i] << "\n";
-                    // std::cout << rhs.view() << "\n";
+                std::cout << "After wake\n";
+                std::cout << rhs.view() << "\n";
+
+                for (i64 m = 0; m < m_assembly.num_surfaces(); m++) {
+                    std::cout << gamma_wake.views()[m].slice(All, Range{-2-i, -1}) << "\n";
                 }
             }
 
@@ -743,6 +759,13 @@ void UVLM_2DOF::run(const UVLM_2DOF_Vars& vars, f32 t_final_nd) {
                     gamma_wing_i.slice(All, -1).to(gamma_wake_i.slice(All, -1-(i+1)));
 
                     // Compute delta
+                    if (m == 1) { // gross hack for missing connectivity on the flap trailing edge
+                        backend->blas->axpy(
+                            -1.0f,
+                            gamma_wing.views()[0].slice(All, -1),
+                            gamma_wing_delta_i.slice(All, 0)
+                        );
+                    }
                     backend->blas->axpy(
                         -1.0f,
                         gamma_wing_i.slice(All, Range{0, -2}),

tests/uvlm_3dof.cpp

@@ -33,8 +33,8 @@
 using namespace vlm;
 using namespace linalg::ostream_overloads;
 
-// #define COUPLED 1
-// #define FREEPLAY_JOINT
+#define COUPLED 1
+#define FREEPLAY_JOINT
 
 // Notes:
 // 1. Nonlinear term is defined twice (structural initialization and in the coupling loop)
@@ -80,9 +80,9 @@ constexpr i32 DOF = 3;
 
 class UVLM_3DOF final: public Simulation {
     public:
-        UVLM_3DOF(const std::string& backend_name, const std::vector<std::string>& meshes);
+        UVLM_3DOF(const std::string& backend_name, const Assembly& assembly);
         ~UVLM_3DOF() = default;
-        void run(const Assembly& assembly, const Vars& vars, f32 t_final);
+        void run(const Vars& vars, f32 t_final);
 
         MultiTensor3D<Location::Device> colloc_d{backend->memory.get()};
         MultiTensor3D<Location::Host> colloc_h{backend->memory.get()};
@@ -134,6 +134,8 @@ class UVLM_3DOF final: public Simulation {
         Tensor1D<Location::Host> dF_h{backend->memory.get()}; // dof
 
         NewmarkBeta integrator{backend.get()};
+
+        Assembly m_assembly;
     private:
         void alloc_buffers();
         void update_wake_and_gamma(i64 iteration);
@@ -144,8 +146,8 @@ class UVLM_3DOF final: public Simulation {
         f32 wing_h(f32 elastic_dst); // get wing elastic height at given dist
 };
 
-UVLM_3DOF::UVLM_3DOF(const std::string& backend_name, const std::vector<std::string>& meshes) : Simulation(backend_name, meshes, false) {
-    ASSERT_EQ(meshes.size(), 2);
+UVLM_3DOF::UVLM_3DOF(const std::string& backend_name, const Assembly& assembly) : Simulation(backend_name, assembly.mesh_filenames(), false), m_assembly(assembly) {
+    ASSERT_EQ(assembly.num_surfaces(), 2);
     solver = backend->create_lu_solver();
     accel_solver = backend->create_lu_solver();
     alloc_buffers();
@@ -244,7 +246,7 @@ void UVLM_3DOF::alloc_buffers() {
 }
 
 #ifdef FREEPLAY_JOINT
-f32 joint_(f32 alpha, f32 M0 = 0.0f, f32 Mf = 0.0f, f32 delta = to_radians(4.24f), f32 a_f = to_radians(-2.12f)) {
+f32 joint(f32 alpha, f32 M0 = 0.0f, f32 Mf = 0.0f, f32 delta = to_radians(4.24f), f32 a_f = to_radians(-2.12f)) {
     if (alpha < a_f) {
         return M0 + alpha - a_f;
     } else if (alpha >= a_f && alpha <= (a_f + delta)) {
@@ -283,6 +285,13 @@ void UVLM_3DOF::update_wake_and_gamma(i64 iteration) {
         }
 
         // Compute delta
+        if (m == 1) { // gross hack for missing connectivity on the flap trailing edge
+            backend->blas->axpy(
+                -1.0f,
+                gamma_wing.views()[0].slice(All, -1),
+                gamma_wing_delta_i.slice(All, 0)
+            );
+        }
         backend->blas->axpy(
             -1.0f,
             gamma_wing_i.slice(All, Range{0, -2}),
@@ -424,7 +433,7 @@ void nondimensionalize_verts(TensorView3D<Location::Host>& verts, const f32 b) {
     }
 }
 
-void UVLM_3DOF::run(const Assembly& assembly, const Vars& v, f32 t_final_nd) {
+void UVLM_3DOF::run(const Vars& v, f32 t_final_nd) {
     const f32 V = v.U / (v.b * v.omega_alpha);
     const f32 mu = v.m / (PI_f * v.rho * v.b*v.b);
 
@@ -456,8 +465,8 @@ void UVLM_3DOF::run(const Assembly& assembly, const Vars& v, f32 t_final_nd) {
     const f32 dy = verts_first_wing(0, -1, 1) - verts_first_wing(0, -2, 1);
     const f32 dz = verts_first_wing(0, -1, 2) - verts_first_wing(0, -2, 2);
     const f32 last_panel_chord = std::sqrt(dx*dx + dy*dy + dz*dz);
-    ASSERT_NEAR(linalg::length(assembly.kinematics()->linear_velocity(0.0f, {0.f, 0.f, 0.f})), v.U / (v.b * v.omega_alpha), 1e-6f);
-    const f32 dt_nd = last_panel_chord / linalg::length(assembly.kinematics()->linear_velocity(0.0f, {0.f, 0.f, 0.f}));
+    ASSERT_NEAR(linalg::length(m_assembly.kinematics()->linear_velocity(0.0f, {0.f, 0.f, 0.f})), v.U / (v.b * v.omega_alpha), 1e-6f);
+    const f32 dt_nd = last_panel_chord / linalg::length(m_assembly.kinematics()->linear_velocity(0.0f, {0.f, 0.f, 0.f}));
     const i64 t_steps = static_cast<i64>(t_final_nd / dt_nd);
     t_h.init({t_steps-2});
 
@@ -489,12 +498,12 @@ void UVLM_3DOF::run(const Assembly& assembly, const Vars& v, f32 t_final_nd) {
     backend->wake_shed(verts_wing.views(), verts_wake.views(), 0);
     // Technically in the initial velocities calculation we should also take into account the IC of the eq of motion
     // In our case, since the IC is an initial position, we can ignore it. 
-    wing_velocities(assembly.surface_kinematics()[0]->transform_dual(0.0f), colloc_h.views()[0], velocities_h.views()[0], velocities.views()[0]);
-    wing_velocities(assembly.surface_kinematics()[0]->transform_dual(0.0f), colloc_h.views()[1], velocities_h.views()[1], velocities.views()[1]);
+    wing_velocities(m_assembly.surface_kinematics()[0]->transform_dual(0.0f), colloc_h.views()[0], velocities_h.views()[0], velocities.views()[0]);
+    wing_velocities(m_assembly.surface_kinematics()[0]->transform_dual(0.0f), colloc_h.views()[1], velocities_h.views()[1], velocities.views()[1]);
 
     rhs.view().fill(0.f);
     backend->rhs_assemble_velocities(rhs.view(), normals_d.views(), velocities.views());
-    backend->rhs_assemble_wake_influence(rhs.view(), gamma_wake.views(), colloc_d.views(), normals_d.views(), verts_wake.views(), 0);
+    backend->rhs_assemble_wake_influence(rhs.view(), gamma_wake.views(), colloc_d.views(), normals_d.views(), verts_wake.views(), m_assembly.lifting(), 0);
     solver->solve(lhs.view(), rhs.view().reshape(rhs.size(), 1));
     update_wake_and_gamma(0);
 
@@ -608,7 +617,7 @@ void UVLM_3DOF::run(const Assembly& assembly, const Vars& v, f32 t_final_nd) {
             // std::printf("iter: %d | h: %f | alpha: %f | h_dot: %f | alpha_dot: %f\n", iteration, u_h.view()(0, i+1), u_h.view()(1, i+1), v_h.view()(0, i+1), v_h.view()(1, i+1));
 
             // Manually compute the transform (update_transforms(assembly, t+dt))
-            auto fs = assembly.kinematics()->transform_dual(t_nd+dt_nd);
+            auto fs = m_assembly.kinematics()->transform_dual(t_nd+dt_nd);
             auto wing_transform_dual = linalg::mul(fs, linalg::mul(heave.transform_dual(0.0f), wing_pitch.transform_dual(0.0f)));
             auto flap_transform_dual = linalg::mul(wing_transform_dual, flap_pitch.transform_dual(0.0f));
             auto wing_transform = dual_to_float(wing_transform_dual);
@@ -631,7 +640,7 @@ void UVLM_3DOF::run(const Assembly& assembly, const Vars& v, f32 t_final_nd) {
 
             rhs.view().fill(0.f);
             backend->rhs_assemble_velocities(rhs.view(), normals_d.views(), velocities.views());
-            backend->rhs_assemble_wake_influence(rhs.view(), gamma_wake.views(), colloc_d.views(), normals_d.views(), verts_wake.views(), i+1);
+            backend->rhs_assemble_wake_influence(rhs.view(), gamma_wake.views(), colloc_d.views(), normals_d.views(), verts_wake.views(), m_assembly.lifting(), i+1);
             solver->solve(lhs.view(), rhs.view());
 
             {
@@ -652,6 +661,13 @@ void UVLM_3DOF::run(const Assembly& assembly, const Vars& v, f32 t_final_nd) {
                     gamma_wing_i.slice(All, -1).to(gamma_wake_i.slice(All, -1-(i+1)));
 
                     // Compute delta
+                    if (m == 1) { // gross hack for missing connectivity on the flap trailing edge
+                        backend->blas->axpy(
+                            -1.0f,
+                            gamma_wing.views()[0].slice(All, -1),
+                            gamma_wing_delta_i.slice(All, 0)
+                        );
+                    }
                     backend->blas->axpy(
                         -1.0f,
                         gamma_wing_i.slice(All, Range{0, -2}),
@@ -662,7 +678,7 @@ void UVLM_3DOF::run(const Assembly& assembly, const Vars& v, f32 t_final_nd) {
                 }
             }
 
-            const linalg::float3 freestream = -assembly.kinematics()->linear_velocity(t_nd+dt_nd, {0.f, 0.f, 0.f});
+            const linalg::float3 freestream = -m_assembly.kinematics()->linear_velocity(t_nd+dt_nd, {0.f, 0.f, 0.f});
             backend->forces_unsteady_multibody(
                 verts_wing.views(),
                 gamma_wing_delta.views(),
@@ -701,12 +717,9 @@ void UVLM_3DOF::run(const Assembly& assembly, const Vars& v, f32 t_final_nd) {
             );
 
             // Aero forces
-            // F_h.view()(0, i+1) = - V*V*cl / (PI_f*mu);
-            // F_h.view()(1, i+1) = 2.f*V*V*cm_a.y / (PI_f*mu);
-            // F_h.view()(2, i+1) = 2.f*V*V*cm_b.y / (PI_f*mu);
-            F_h.view()(0, i+1) = - 5.f * cl / (PI_f*mu);
+            F_h.view()(0, i+1) = - V*V * cl / (PI_f*mu);
             F_h.view()(1, i+1) = 2.f*V*V*cm_a.y / (PI_f*mu);
-            F_h.view()(2, i+1) = 2.f*V*V*cm_b.y / (PI_f*mu);
+            F_h.view()(2, i+1) = pow(1-v.c, 2)*V*V*cm_b.y / (2.0f * PI_f*mu);
 
             // deltaF for rhs
             const f32 sigma = v.omega_h / v.omega_alpha;
@@ -734,10 +747,11 @@ void UVLM_3DOF::run(const Assembly& assembly, const Vars& v, f32 t_final_nd) {
 
 int main() {
     // vlm::Executor::instance(1);
-    const std::vector<std::vector<std::string>> meshes = {{
-        "../../../../mesh/3dof_wing_9x5.x",
-        "../../../../mesh/3dof_flap_3x5.x"
-    }}; // vector of vector meh
+    std::vector<std::vector<std::pair<std::string, bool>>> meshes;
+    meshes.push_back({
+        {"../../../../mesh/3dof_wing_9x5.x", false},
+        {"../../../../mesh/3dof_flap_3x5.x", true}
+    });
     const std::vector<std::string> backends = {"cpu"};
 
     auto simulations = tiny::make_combination(meshes, backends);
@@ -772,7 +786,7 @@ int main() {
     v.U = 12.36842f; // m/s
     // v.U = 5.0f; // m/s
 
-    const f32 t_final_nd = 3.f * v.omega_alpha;
+    const f32 t_final_nd = 10.f * v.omega_alpha;
 
     KinematicsTree kinematics_tree;
 
@@ -784,11 +798,13 @@ int main() {
         });
     });
 
-    for (const auto& [mesh_names, backend_name] : simulations) {
+    for (const auto& [surface, backend_name] : simulations) {
         Assembly assembly(freestream);
-        assembly.add(mesh_names.get()[0], freestream);
-        UVLM_3DOF simulation{backend_name, mesh_names};
-        simulation.run(assembly, v, t_final_nd);
+        for (const auto& [mesh_name, lifting] : surface.get()) {
+            assembly.add(mesh_name, freestream, lifting);
+        }
+        UVLM_3DOF simulation{backend_name, assembly};
+        simulation.run(v, t_final_nd);
     }
 
     return 0;

vlm/backends/cpu/src/vlm_backend_cpu.cpp

@@ -322,39 +322,45 @@ void BackendCPU::rhs_assemble_velocities(TensorView<f32, 1, Location::Device>& r
 
 void BackendCPU::rhs_assemble_wake_influence(TensorView<f32, 1, Location::Device>& rhs, const MultiTensorView2D<Location::Device>& gamma_wake, const MultiTensorView3D<Location::Device>& colloc, const MultiTensorView3D<Location::Device>& normals, const MultiTensorView3D<Location::Device>& verts_wake, const std::vector<bool>& lifting, i32 iteration) {
     // const tiny::ScopedTimer timer("Wake Influence");
+    assert(lifting.size() == normals.size());
 
     tf::Taskflow taskflow;
 
     auto begin = taskflow.placeholder();
     auto end = taskflow.placeholder();
 
-    auto wake_influence = taskflow.for_each_index((i64)0, rhs.size(), [&] (i64 idx) {
-        // Loop over the wakes
-        for (i32 i = 0; i < normals.size(); i++) {
-            if (!lifting[i]) continue;
-            const auto& normals_i = normals[i];
-            const auto& colloc_i = colloc[i];
-            const auto& gamma_wake_i = gamma_wake[i];
-            const auto& verts_wake_i = verts_wake[i];
-            ispc::kernel_wake_influence(
-                colloc_i.ptr() + idx,
-                colloc_i.stride(2),
-                normals_i.ptr() + idx,
-                normals_i.stride(2),
-                verts_wake_i.ptr(),
-                verts_wake_i.stride(2),
-                verts_wake_i.shape(1),
-                verts_wake_i.shape(0),
-                gamma_wake_i.ptr(), 
-                rhs.ptr() + idx,
-                sigma_vatistas,
-                iteration
-            );
-        }
-    });
+    i64 offset = 0;
+    for (i32 m_i = 0; m_i < normals.size(); m_i++) {
+        const i64 m_i_num_panels = normals[m_i].shape(0) * normals[m_i].shape(1);
+        auto wake_influence = taskflow.for_each_index((i64)0, m_i_num_panels, [=, &lifting] (i64 idx) {
+            // Loop over the wakes
+            for (i32 m_j = 0; m_j < normals.size(); m_j++) {
+                if (!lifting[m_j]) continue;
+                const auto& normals_i = normals[m_i];
+                const auto& colloc_i = colloc[m_i];
+                const auto& gamma_wake_j = gamma_wake[m_j];
+                const auto& verts_wake_j = verts_wake[m_j];
+                ispc::kernel_wake_influence(
+                    colloc_i.ptr() + idx, // technically incorrect as it doesnt consider the edge case of non-contiguous buffer
+                    colloc_i.stride(2),
+                    normals_i.ptr() + idx, // identical as colloc_i
+                    normals_i.stride(2),
+                    verts_wake_j.ptr(),
+                    verts_wake_j.stride(2),
+                    verts_wake_j.shape(1),
+                    verts_wake_j.shape(0),
+                    gamma_wake_j.ptr(), 
+                    rhs.ptr() + offset + idx,
+                    sigma_vatistas,
+                    iteration
+                );
+            }
+        });
 
-    begin.precede(wake_influence);
-    wake_influence.precede(end);
+        begin.precede(wake_influence);
+        wake_influence.precede(end);
+        offset += m_i_num_panels;
+    }
 
     Executor::get().run(taskflow).wait();
 }

vlm/include/vlm_memory.hpp

@@ -253,7 +253,7 @@ class TensorView {
 namespace detail {
     template<typename T, i64 Dim, Location L>
     void print_tensor_recursive(std::ostream& os, const TensorView<T, Dim, L>& view, 
-                              std::array<i64, Dim>& indices, i64 current_dim = Dim - 1) {
+                              std::array<i64, static_cast<size_t>(Dim)>& indices, i64 current_dim) {
         if (current_dim < 0) {
             os << "(";
             for (i64 i = 0; i < Dim; ++i) {
@@ -273,8 +273,8 @@ namespace detail {
 
 template<typename T, i64 Dim, Location L>
 std::ostream& operator<<(std::ostream& os, const TensorView<T, Dim, L>& view) {
-    std::array<i64, Dim> indices{};
-    detail::print_tensor_recursive(os, view, indices);
+    std::array<i64, static_cast<size_t>(Dim)> indices{};
+    detail::print_tensor_recursive(os, view, indices, Dim - 1);
     return os;
 }