add true vortex panels

TODO

@@ -3,4 +3,5 @@
 - [X] Create additional SolverData struct for vatistas parameter (ideally find another paramater too) -> removed and directly put it as member of backend base class
 - [ ] Rename many functions and harmonize the API
 - [ ] Decide what to do with SoA3D and evaluate the fusion of the 3 buffers into one 
-- [ ] using a view for multidimensional indexing.
+- [ ] using a view for multidimensional indexing.
+- [ ] Instead of view, use a function that returns a tuple of the panel's vertex indices that takes in arg i, j values for chordwise and spanwise location

data/theodorsen.py

@@ -38,7 +38,7 @@ def theo_fun(k):
 a = ar / c # full span
 pitch_axis = -0.5 # leading edge
 
-cycles = 4 # number of periods
+cycles = 3 # number of periods
 
 uvlm_cl = []
 uvlm_t = []
@@ -82,7 +82,7 @@ def theo_fun(k):
     # def heave(t): return 0
 
     # pure heaving
-    def pitch(t): return np.radians(-5)
+    def pitch(t): return 0
     def heave(t): return -amplitude * np.sin(omega * t)
 
     # pure pitching
@@ -118,15 +118,15 @@ def cl_theodorsen(t: float): # using Wagner functions and Kholodar formulation
 error = np.sqrt(np.dot(difference, difference) / (n-uvlm_cycle_idx)) 
 print(f"Error: {100 * error / np.max(np.abs(analytical_cl)):.3f}%", )
 
-katz_t = []
-katz_cl = []
-with open("data/katz/katz_cl.txt", "r") as f:
-    for line in f:
-        t, cl = line.split()
-        katz_t.append(float(t))
-        katz_cl.append(float(cl))
+# katz_t = []
+# katz_cl = []
+# with open("data/katz/katz_cl.txt", "r") as f:
+#     for line in f:
+#         t, cl = line.split()
+#         katz_t.append(float(t))
+#         katz_cl.append(float(cl))
 
-axs["time"].scatter(katz_t[1:], katz_cl[1:], label="Katz", facecolors='none', edgecolors='g', s=20)
+# axs["time"].scatter(katz_t[1:], katz_cl[1:], label="Katz", facecolors='none', edgecolors='g', s=20)
 
 axs["time"].scatter(uvlm_t, uvlm_cl, label="UVLM", facecolors='none', edgecolors='r', s=20)
 # axs["time"].scatter(uvlm_t, uvlm_z, label="UVLM z", facecolors='none', edgecolors='r', s=20)

tests/test_uvlm_theodorsen.cpp

@@ -82,8 +82,8 @@ int main() {
     const u64 ni = 10;
     const u64 nj = 5;
     // vlm::Executor::instance(1);
-    const std::vector<std::string> meshes = {"../../../../mesh/rectangular_5x10.x"};
-    //const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_" + std::to_string(ni) + "x" + std::to_string(nj) + ".x"};
+    //const std::vector<std::string> meshes = {"../../../../mesh/rectangular_5x10.x"};
+    const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_" + std::to_string(ni) + "x" + std::to_string(nj) + ".x"};
     const std::vector<std::string> backends = get_available_backends();
 
     auto solvers = tiny::make_combination(meshes, backends);
@@ -92,24 +92,24 @@ int main() {
     const f32 b = 0.5f; // half chord
 
     // Define simulation length
-    const f32 cycles = 4.0f;
+    const f32 cycles = 3.0f;
     const f32 u_inf = 1.0f; // freestream velocity
     const f32 amplitude = 0.1f; // amplitude of the wing motion
-    const f32 k = 0.5; // reduced frequency
+    const f32 k = 0.75; // reduced frequency
     const f32 omega = k * 2.0f * u_inf / (2*b);
     const f32 t_final = cycles * 2.0f * PI_f / omega; // 4 periods
     //const f32 t_final = 5.0f;
 
     Kinematics kinematics{};
 
-    tmatrix initial_pose = linalg::rotation_matrix(
+    const f32 initial_angle = 0.0f;
+
+    const tmatrix initial_pose = linalg::rotation_matrix(
         linalg::alias::float3{0.0f, 0.0f, 0.0f}, // take into account quarter chord panel offset
         linalg::alias::float3{0.0f, 1.0f, 0.0f},
-        to_radians(-5.f)
+        to_radians(initial_angle)
     );
 
-    //tmatrix initial_pose = linalg::identity;
-
     // Periodic heaving
     kinematics.add([=](f32 t) {
         return linalg::translation_matrix(linalg::alias::float3{
@@ -182,7 +182,6 @@ int main() {
                 for (u64 i = 1; i < trailing_vertices.size; i++) {
                     dt = std::min(dt, segment_chord / kinematics.velocity_magnitude(t, {trailing_vertices.x[i], trailing_vertices.y[i], trailing_vertices.z[i], 1.0f}));
                 }
-                dt = 0.2f;
 
                 auto transform = kinematics.relative_displacement(t, t+dt);
                 for (u64 i = 0; i < trailing_vertices.size; i++) {
@@ -261,7 +260,7 @@ int main() {
                 if (idx == 0) {
                     linalg::alias::float4 freestream_vel = -kinematics.velocity(t, colloc_pt, 1);
                     freestream = {freestream_vel.x, freestream_vel.y, freestream_vel.z};
-                    std::cout << freestream << "\n";
+                    std::cout << "Freestream: " << freestream << "\n";
                     const f32 analytical_vel = - amplitude * omega * std::cos(omega * t);
                     const f32 rel_error = 100.0f * std::abs((analytical_vel - local_velocity.z) / analytical_vel);
                     std::cout << "vel error:" << rel_error << "%\n";

tests/test_vlm_elliptic_coeffs.cpp

@@ -133,7 +133,7 @@ int main(int argc, char **argv) {
             std::printf(">>> Alpha: %.1f | CL = %.6f CD = %.6f CMx = %.6f CMy = %.6f CMz = %.6f\n", to_degrees(flow.alpha), coeffs.cl, coeffs.cd, coeffs.cm.x, coeffs.cm.y, coeffs.cm.z);
             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;
+            if (cl_rerr > 0.03f || cd_rerr > 0.03f) return 1;
         } 
     }
     return 0;

vlm/backends/cpu/src/vlm_backend_cpu.cpp

@@ -282,7 +282,7 @@ f32 BackendCPU::compute_coefficient_unsteady_cl(const linalg::alias::float3& fre
             const f32 gamma_dt = (gamma[li] - gamma_prev[li]) / dt; // backward difference
 
             // Joukowski method
-            force += rho * delta_gamma[li] * linalg::cross(freestream, v1 - v0);
+            force += rho * delta_gamma[li] * linalg::cross(V, v1 - v0);
             force += rho * gamma_dt * mesh.area[li] * normal;
 
             // Katz Plotkin method

vlm/src/vlm_mesh.cpp

@@ -42,6 +42,27 @@ std::unique_ptr<Mesh> vlm::create_mesh(const std::string& filename, const u64 nw
 }
 
 void Mesh::create_vortex_panels() {
+    for (u64 i = 0; i < nc; i++) {
+        for (u64 j = 0; j < ns+1; j++) {
+            const u64 v0 = (i+0) * (ns+1) + j;
+            const u64 v3 = (i+1) * (ns+1) + j;
+
+            v.x[v0] = 0.75f * v.x[v0] + 0.25f * v.x[v3];
+            v.y[v0] = 0.75f * v.y[v0] + 0.25f * v.y[v3];
+            v.z[v0] = 0.75f * v.z[v0] + 0.25f * v.z[v3];
+        }
+    }
+
+    // Trailing edge vertices
+    const u64 i = nc;
+    for (u64 j = 0; j < ns+1; j++) {
+        const u64 v0 = (i+0) * (ns+1) + j;
+        const u64 v3 = (i+1) * (ns+1) + j;
+
+        v.x[v3] = 1.25f * v.x[v3] - 0.25f * v.x[v0];
+        v.y[v3] = 1.25f * v.y[v3] - 0.25f * v.y[v0];
+        v.z[v3] = 1.25f * v.z[v3] - 0.25f * v.z[v0];
+    }
 }
 
 void Mesh::init() {