added core uvlm algorithm

good progress

TODO

@@ -0,0 +1,4 @@
+- Get rid of 3D vectors and use homogeneous coordinates everywhere (4d vector).
+- Move some of the logic from theodorsen test case to its own class
+- Create additional SolverData struct for vatistas parameter (ideally find another paramater too)
+- Rename many functions and harmonize the API

headeronly/tinytypes.hpp

@@ -1,5 +1,7 @@
 #include <cstdint>
 #include <cmath>
+#include <cstring>
+#include <vector>
 
 namespace tiny {
 

tests/test_uvlm_theodorsen.cpp

@@ -8,6 +8,7 @@
 #include "tinycombination.hpp"
 
 #include "vlm.hpp"
+#include "vlm_data.hpp"
 #include "vlm_utils.hpp"
 
 using namespace vlm;
@@ -37,13 +38,17 @@ class Kinematics {
         return linalg::mul(displacement(t1), linalg::inverse(displacement(t0)));
     }
 
-    f32 absolute_velocity(f32 t, const linalg::alias::float4& vertex) {
-        return linalg::length((linalg::mul(relative_displacement(t, t+EPS_sqrt_f), vertex)-vertex)/EPS_sqrt_f);
+    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;
+    }
+
+    f32 velocity_magnitude(f32 t, const linalg::alias::float4& vertex) {
+        return linalg::length(velocity(t, vertex));
     }
 };
 
 int main() {
-    const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_5x200.x"};
+    const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_2x8.x"};
     const std::vector<std::string> backends = get_available_backends();
 
     auto solvers = tiny::make_combination(meshes, backends);
@@ -71,29 +76,39 @@ int main() {
     });
 
     for (const auto& [mesh_name, backend_name] : solvers) {
-        UVLM solver{};
-        solver.mesh = create_mesh(mesh_name);
+        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 * (solver.mesh->chord_length(0) + solver.mesh->chord_length(1))) / kinematics.absolute_velocity(t, {0.f, 0.f, 0.f, 1.f});
+            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});
             t += dt;
             vec_t.push_back(t);
         }
 
-        solver.mesh->resize_wake(vec_t.size()-1); // +1 for the initial pose
-        solver.backend = create_backend(backend_name, *solver.mesh); // create after mesh has been resized
-
-        // f32 t = 0.0f;
-        // for (u64 i = 0; i < vec_dt.size(); i++) {
-        //     FlowData flow{};
-        //     solver.solve(flow);
-        //     solver.mesh->shed_wake();
-        //     solver.mesh->move(linalg::mul(displacement(t+vec_dt[i]), linalg::inverse(displacement(t))));
-        //     t += vec_dt[i];
-        // }
+        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
+
+        // Precompute the LHS since wing geometry is constant
+        FlowData flow_dummy{0.0f, 0.0f, 1.0f, 1.0f};
+        backend->compute_lhs(flow_dummy);
+        backend->lu_factor();
+        for (u64 i = 0; i < vec_t.size()-1; i++) {
+            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};
+            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]));
+            backend->shed_gamma();
+        }
     }
 
     return 0;

tests/test_vlm_elliptic_coeffs.cpp

@@ -134,8 +134,7 @@ int main(int argc, char **argv) {
             std::printf(">>> Analytical CL: %.6f | Abs Error: %.3E | Relative Error: %.5f%% \n", analytical_cl, cl_aerr, cl_rerr*100.f);
             std::printf(">>> Analytical CD: %.6f | Abs Error: %.3E | Relative Error: %.5f%% \n", analytical_cd, cd_aerr, cd_rerr*100.f);
             if (cl_rerr > 0.03f || cd_rerr > 0.01f) return 1;
-        }
+        } 
     }
-
     return 0;
 }

vlm/backends/cpu/include/vlm_backend_cpu.hpp

@@ -8,9 +8,10 @@ namespace vlm {
 class BackendCPU : public Backend {
     public:
         std::vector<f32> lhs;
+        std::vector<f32> wake_buffer; // (ns*nc) * ns
         std::vector<f32> rhs;
         std::vector<i32> ipiv;
-        std::vector<f32> gamma;
+        std::vector<f32> gamma; // ncw * ns
         std::vector<f32> delta_gamma;
         std::vector<f32> trefftz_buffer;
 
@@ -20,6 +21,8 @@ class BackendCPU : public Backend {
         void compute_lhs(const FlowData& flow) override;
         void compute_rhs(const FlowData& flow) override;
         void compute_rhs(const FlowData& flow, const std::vector<f32>& section_alphas) override; 
+        void add_wake_influence(const FlowData& flow) override;
+        void shed_gamma() override;
         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;

vlm/backends/cpu/src/vlm_backend_cpu.cpp

@@ -22,16 +22,19 @@ using namespace vlm;
 
 BackendCPU::~BackendCPU() = default; // Destructor definition
 
+// TODO: replace any kind of size arithmetic with methods
 BackendCPU::BackendCPU(Mesh& mesh) : Backend(mesh) {
     lhs.resize(mesh.nb_panels_wing() * mesh.nb_panels_wing());
+    wake_buffer.resize(mesh.nb_panels_wing() * mesh.ns);
     rhs.resize(mesh.nb_panels_wing());
     ipiv.resize(mesh.nb_panels_wing());
-    gamma.resize(mesh.nb_panels_wing());
+    gamma.resize((mesh.nc + mesh.nw) * mesh.ns); // store wake gamma as well
     delta_gamma.resize(mesh.nb_panels_wing());
     trefftz_buffer.resize(mesh.ns);
 }
 
 void BackendCPU::reset() {
+    // TODO: verify that these are needed
     std::fill(gamma.begin(), gamma.end(), 0.0f);
     std::fill(lhs.begin(), lhs.end(), 0.0f);
     std::fill(rhs.begin(), rhs.end(), 0.0f);
@@ -71,7 +74,7 @@ void BackendCPU::compute_lhs(const FlowData& flow) {
 
     u64 idx = m.nc - 1;
     auto cond = taskflow.emplace([&]{
-        return idx < m.nc + m.nw ? 0 : 1; // 0 means continue, 1 means break
+        return idx < m.nc + m.current_nw ? 0 : 1; // 0 means continue, 1 means break
     });
     auto wake_pass = taskflow.for_each_index(zero, m.ns, [&] (u64 j) {
         const u64 ia = (m.nc - 1) * m.ns + j;
@@ -106,6 +109,42 @@ void BackendCPU::compute_rhs(const FlowData& flow) {
     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");
+
+    tf::Taskflow taskflow;
+
+    Mesh& m = mesh;
+    ispc::MeshProxy mesh_proxy = {
+        m.ns, m.nc, m.nb_panels_wing(),
+        {m.v.x.data(), m.v.y.data(), m.v.z.data()}, 
+        {m.colloc.x.data(), m.colloc.y.data(), m.colloc.z.data()},
+        {m.normal.x.data(), m.normal.y.data(), m.normal.z.data()}
+    };
+
+    // 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);
+        }
+
+        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);
+    }
+}
+
+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() + m.ns * (m.nc-1), gamma.data() + m.nb_panels_wing(), gamma.data() + wake_row_start);
+}
+
 void BackendCPU::compute_rhs(const FlowData& flow, const std::vector<f32>& section_alphas) {
     tiny::ScopedTimer timer("Rebuild RHS");
     assert(section_alphas.size() == mesh.ns);

vlm/backends/cuda/src/vlm_backend_cuda.cu

@@ -312,7 +312,7 @@ void BackendCUDA::compute_lhs(const FlowData& flow) {
     // CHECK_CUDA(cudaDeviceSynchronize());
 
     dim3 grid_size2(get_grid_size(mesh.nb_panels_wing(), block_size.x), get_grid_size(mesh.ns, block_size.y));
-    for (u64 offset = 0; offset < mesh.nw + 1; offset++) {
+    for (u64 offset = 0; offset < mesh.current_nw + 1; offset++) {
         kernel_influence_cuda<<<grid_size2, block_size>>>(d_mesh, d_lhs, (mesh.nc - 1) * mesh.ns, mesh.ns, offset*mesh.ns, flow.sigma_vatistas);
         // cudaError_t error = cudaGetLastError();
         // if (error != cudaSuccess) {

vlm/include/vlm_backend.hpp

@@ -14,6 +14,8 @@ class Backend {
         virtual void compute_lhs(const FlowData& flow) = 0;
         virtual void compute_rhs(const FlowData& flow) = 0;
         virtual void compute_rhs(const FlowData& flow, const std::vector<f32>& section_alphas) = 0; 
+        virtual void add_wake_influence(const FlowData& flow) = 0;
+        virtual void shed_gamma() = 0;
         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;

vlm/include/vlm_data.hpp

@@ -9,6 +9,7 @@ namespace vlm {
 // Flow characteristics
 class FlowData {
     public:
+    // TODO: maybe move these to methods
     const f32 alpha; // global angle of attack
     const f32 beta; // global sideslip angle
     const f32 u_inf; // freestream velocity magnitude
@@ -20,12 +21,13 @@ class FlowData {
     const linalg::alias::float3 stream_axis;
 
     FlowData(const f32 alpha_, const f32 beta_, const f32 u_inf_, const f32 rho_);
+    FlowData(const linalg::alias::float3& freestream_, const f32 rho_);
 };
 
-class PoseData {
+class SolverData {
     public:
     // Poses defined in global stationary frame
-    const linalg::alias::float4x4 body; // plane body pose
+    const f32 sigma_vatistas = 0.0f; // vatistas coefficient
 };
 
 class AeroCoefficients {

vlm/src/vlm_data.cpp

@@ -10,4 +10,14 @@ FlowData::FlowData(const f32 alpha_, const f32 beta_, const f32 u_inf_, const f3
     freestream(compute_freestream(u_inf, alpha, beta)),
     lift_axis(compute_lift_axis(freestream)),
     stream_axis(compute_stream_axis(freestream)) {
+}
+
+FlowData::FlowData(const linalg::alias::float3& freestream_, const f32 rho_): 
+    alpha(std::atan(freestream_.z / std::sqrt(freestream_.x*freestream_.x + freestream_.y*freestream_.y))),
+    beta(std::atan(freestream_.y / freestream_.x)),
+    u_inf(linalg::length(freestream_)),
+    rho(rho_),
+    freestream(freestream_),
+    lift_axis(compute_lift_axis(freestream)),
+    stream_axis(compute_stream_axis(freestream)) {
 }

vlm/src/vlm_mesh.cpp

@@ -206,6 +206,7 @@ void Mesh::update_wake(const linalg::alias::float3& freestream) {
         v.z[i + v_ns] = v.z[i] + off_z;
     }
     compute_metrics_wake();
+    current_nw = 1;
 }
 
 // https://publications.polymtl.ca/2555/1/2017_MatthieuParenteau.pdf (Eq 3.4 p21)
@@ -392,6 +393,8 @@ void Mesh::move(const linalg::alias::float4x4& transform) {
     std::copy(v.x.data() + src_begin, v.x.data() + src_end, v.x.data() + dst_begin - (ns + 1));
     std::copy(v.y.data() + src_begin, v.y.data() + src_end, v.y.data() + dst_begin - (ns + 1));
     std::copy(v.z.data() + src_begin, v.z.data() + src_end, v.z.data() + dst_begin - (ns + 1));
+
+    current_nw++;
 }
 
 void Mesh::resize_wake(const u64 nw_) {