added rollup

single threaded for the moment. It's super cool !

data/displacement_real.py

@@ -20,7 +20,7 @@ def coords_to_verts(coords, nc, ns):
 coords_wing = []
 coords_wake = []
 
-with open("build/windows/x64/debug/wing_data.txt") as f:
+with open("build/windows/x64/release/wing_data.txt") as f:
     # first line
     nc, ns = map(int, f.readline().split())
     timesteps = int(f.readline())
@@ -32,7 +32,7 @@ def coords_to_verts(coords, nc, ns):
         coords[2, :] = np.array(list(map(float, f.readline().split())))
         coords_wing.append(coords)
 
-with open("build/windows/x64/debug/wake_data.txt") as f:
+with open("build/windows/x64/release/wake_data.txt") as f:
     # first line
     f.readline() # skip first line
     f.readline() # skip second line

data/theodorsen.py

@@ -65,12 +65,12 @@ def atime(t: float): return 2. * u_inf * t / c
     # def heave(t): return 0
 
     # pure heaving
-    # def pitch(t): return 0
-    # def heave(t): return -amplitude * np.sin(omega * t)
+    def pitch(t): return 0
+    def heave(t): return -amplitude * np.sin(omega * t)
 
     # pure pitching
-    def pitch(t): return np.radians(np.sin(omega * t))
-    def heave(t): return 0
+    # def pitch(t): return np.radians(np.sin(omega * t))
+    # def heave(t): return 0
 
     def w(s: float): 
         return u_inf * pitch(s) + derivative(heave, s) + b * (0.5 - pitch_axis) * derivative(pitch, s)
@@ -100,7 +100,7 @@ def cl_theodorsen(t: float): # using Wagner functions and Kholodar formulation
 uvlm_t = []
 uvlm_z = []
 uvlm_alpha = []
-with open("build/windows/x64/debug/cl_data.txt", "r") as f:
+with open("build/windows/x64/release/cl_data.txt", "r") as f:
     for line in f:
         time_step, z, cl, alpha = line.split()
         uvlm_t.append(float(time_step))

tests/test_uvlm_theodorsen.cpp

@@ -77,7 +77,7 @@ int main() {
     const tiny::ScopedTimer timer("UVLM TOTAL");
     // vlm::Executor::instance(1);
     // const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_2x8.x"};
-    const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_2x2.x"};
+    const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_5x10.x"};
 
     const std::vector<std::string> backends = get_available_backends();
 
@@ -239,6 +239,7 @@ int main() {
                 cl_data << t << " " << mesh->v.z[0] << " " << cl_unsteady << " " << std::sin(omega * t) << "\n";
                 #endif
             }
+            backend->wake_rollup(dt);
             backend->shed_gamma(); // shed before moving & incrementing currentnw
             mesh->move(kinematics.relative_displacement(t, t+dt));
         }

vlm/backends/cpu/include/vlm_backend_cpu.hpp

@@ -29,6 +29,7 @@ class BackendCPU : public Backend {
         void compute_rhs(const FlowData& flow) override;
         void compute_rhs(const SoA_3D_t<f32>& velocities) override; 
         void add_wake_influence() override;
+        void wake_rollup(float dt) override;
         void shed_gamma() override;
         void lu_factor() override;
         void lu_solve() override;

vlm/backends/cpu/src/vlm_backend_cpu.cpp

@@ -201,6 +201,25 @@ void BackendCPU::add_wake_influence() {
     // Executor::get().run(taskflow).wait();
 }
 
+void BackendCPU::wake_rollup(float dt) {
+    const tiny::ScopedTimer timer("Wake Rollup");
+
+    ispc::MeshView mesh_view = {
+        mesh.nc, mesh.ns, mesh.nw, mesh.current_nw,
+        {mesh.v.x.data(), mesh.v.y.data(), mesh.v.z.data()},
+        {mesh.colloc.x.data(), mesh.colloc.y.data(), mesh.colloc.z.data()},
+        {mesh.normal.x.data(), mesh.normal.y.data(), mesh.normal.z.data()},
+    };
+
+    const u64 wake_vertices_begin = (mesh.nc + mesh.nw - mesh.current_nw + 1) * (mesh.ns+1);
+    const u64 wake_vertices_end = (mesh.nc + mesh.nw + 1) * (mesh.ns + 1);
+
+    // parallel for
+    for (u64 vidx = wake_vertices_begin; vidx < wake_vertices_end; vidx++) {
+        ispc::kernel_induced_vel(mesh_view, vidx, dt, gamma.data(), sigma_vatistas);
+    }
+}
+
 void BackendCPU::shed_gamma() {
     const Mesh& m = mesh;
     const u64 wake_row_start = (m.nc + m.nw - m.current_nw - 1) * m.ns;

vlm/backends/cpu/src/vlm_backend_cpu_kernels.ispc

@@ -43,14 +43,25 @@ export struct MeshProxy {
     uniform SoA3D normal; // normals
 };
 
+export struct MeshView {
+    uniform uint64 ni; // wing chord panels
+    uniform uint64 nj; // wing span panels
+    uniform uint64 nw; // wake chord capacity
+    uniform uint64 cw; // number of active wake panels
+    uniform SoA3D v; // vertices
+    uniform SoA3D colloc; // collocation points
+    uniform SoA3D normal; // normals
+};
+
 // Bio-savart Kernel
 #define RCUT 1e-10f
 #define RCUT2 1e-5f
 
 #define PI_f 3.141593f
 
-inline float3 kernel_biosavart(float3& colloc, const uniform float3& vertex1, const uniform float3& vertex2, const uniform float& sigma) {
-    uniform float3 r0 = vertex2 - vertex1;
+template<typename C, typename V>
+inline float3 kernel_biosavart(C& colloc, const V& vertex1, const V& vertex2, const uniform float& sigma) {
+    V r0 = vertex2 - vertex1;
     float3 r1 = colloc - vertex1;
     float3 r2 = colloc - vertex2;
     // Katz Plotkin, Low speed Aero | Eq 10.115
@@ -64,13 +75,14 @@ inline float3 kernel_biosavart(float3& colloc, const uniform float3& vertex1, co
         return res;
     }
 
-    uniform float smoother = sigma*sigma*length2(r0);
+    float smoother = sigma*sigma*length2(r0);
 
     float coeff = (dot(r0,r1)*r2_norm - dot(r0, r2)*r1_norm) / (4.0f*PI_f*sqrt(square*square + smoother*smoother)*r1_norm*r2_norm);
     return r1r2cross * coeff;
 }
 
-inline void kernel_symmetry(float3& inf, float3 colloc, const uniform float3& vertex0, const uniform float3& vertex1, const uniform float& sigma) {
+template<typename C, typename V>
+inline void kernel_symmetry(float3& inf, C colloc, const V& vertex0, const V& vertex1, const uniform float& sigma) {
     float3 induced_speed = kernel_biosavart(colloc, vertex0, vertex1, sigma);
     inf.x += induced_speed.x;
     inf.y += induced_speed.y;
@@ -135,7 +147,6 @@ export void kernel_influence2(
     const uniform uint64 v3 = v0 + m.ns+1;
     const uniform uint64 v2 = v3 + 1;
 
-
     // Broadcast vertices
     uniform float3 vertex0 = {m.v.x[v0], m.v.y[v0], m.v.z[v0]};
     uniform float3 vertex1 = {m.v.x[v1], m.v.y[v1], m.v.z[v1]};
@@ -169,6 +180,68 @@ export void kernel_influence2(
     }
 }
 
+export void kernel_induced_vel(
+    uniform const MeshView& m, uniform uint64 vidx, uniform float dt,
+    uniform float* uniform gamma, uniform float sigma
+    ) {
+
+    const uniform float3 vertex_influenced = {m.v.x[vidx], m.v.y[vidx], m.v.z[vidx]};
+
+    float3 induced_vel = {0.0f, 0.0f, 0.0f};
+
+    // Wing influence
+    for (uint64 i = 0; i < m.ni; i++) {
+        foreach(j = 0 ... m.nj) {
+            uint64 v0 = (i+0) * (m.nj+1) + j;
+            uint64 v1 = (i+0) * (m.nj+1) + j + 1;
+            uint64 v2 = (i+1) * (m.nj+1) + j + 1;
+            uint64 v3 = (i+1) * (m.nj+1) + j;
+
+            // Loads
+            const float3 vertex0 = {m.v.x[v0], m.v.y[v0], m.v.z[v0]};
+            const float3 vertex1 = {m.v.x[v1], m.v.y[v1], m.v.z[v1]};
+            const float3 vertex2 = {m.v.x[v2], m.v.y[v2], m.v.z[v2]};
+            const float3 vertex3 = {m.v.x[v3], m.v.y[v3], m.v.z[v3]};
+
+            float3 ind = {0.0f, 0.0f, 0.0f};
+            kernel_symmetry(ind, vertex_influenced, vertex0, vertex1, sigma);
+            kernel_symmetry(ind, vertex_influenced, vertex1, vertex2, sigma);
+            kernel_symmetry(ind, vertex_influenced, vertex2, vertex3, sigma);
+            kernel_symmetry(ind, vertex_influenced, vertex3, vertex0, sigma);
+
+            induced_vel += ind * gamma[i * m.nj + j];
+        }
+    }
+
+    // Wake influence
+    for (uint64 i = m.ni + m.nw - m.cw; i < m.ni + m.nw; i++) {
+        foreach(j = 0 ... m.nj) {
+            uint64 v0 = (i+0) * (m.nj+1) + j;
+            uint64 v1 = (i+0) * (m.nj+1) + j + 1;
+            uint64 v2 = (i+1) * (m.nj+1) + j + 1;
+            uint64 v3 = (i+1) * (m.nj+1) + j;
+
+            // Loads
+            const float3 vertex0 = {m.v.x[v0], m.v.y[v0], m.v.z[v0]};
+            const float3 vertex1 = {m.v.x[v1], m.v.y[v1], m.v.z[v1]};
+            const float3 vertex2 = {m.v.x[v2], m.v.y[v2], m.v.z[v2]};
+            const float3 vertex3 = {m.v.x[v3], m.v.y[v3], m.v.z[v3]};
+
+            float3 ind = {0.0f, 0.0f, 0.0f};
+            kernel_symmetry(ind, vertex_influenced, vertex0, vertex1, sigma);
+            kernel_symmetry(ind, vertex_influenced, vertex1, vertex2, sigma);
+            kernel_symmetry(ind, vertex_influenced, vertex2, vertex3, sigma);
+            kernel_symmetry(ind, vertex_influenced, vertex3, vertex0, sigma);
+
+            induced_vel += ind * gamma[i * m.nj + j];
+        }
+    }
+
+    m.v.x[vidx] += reduce_add(induced_vel.x) * dt;
+    m.v.y[vidx] += reduce_add(induced_vel.y) * dt;
+    m.v.z[vidx] += reduce_add(induced_vel.z) * dt;
+}
+
 export uniform float kernel_trefftz_cd(
     uniform const MeshProxy& m,
     uniform float* uniform gamma,

vlm/include/vlm_backend.hpp

@@ -16,6 +16,7 @@ class Backend {
         virtual void compute_rhs(const FlowData& flow) = 0;
         virtual void compute_rhs(const SoA_3D_t<f32>& velocities) = 0; 
         virtual void add_wake_influence() = 0;
+        virtual void wake_rollup(float dt) = 0;
         virtual void shed_gamma() = 0;
         virtual void lu_factor() = 0;
         virtual void lu_solve() = 0;