experimental true quarter chord panels

doesnt offer any gains in accuracy tho

config/elliptic.vlm

@@ -8,7 +8,7 @@ alphas = [1, 3, 5]
 s_ref = 3.92699 # half wing area
 c_ref = 0.85
 sigma_vatistas = 0.0
-backend = cuda
+backend = cpu
 
 <mesh>
 wake_included = false

tests/test_linear.cpp

@@ -112,7 +112,7 @@ int main(int argc, char **argv) {
 
     const std::vector<std::string> meshes = {"../../../../mesh/elliptic_64x64.x"};
     const std::vector<std::string> backends = get_available_backends();
-    std::vector<f32> test_alphas = {0, 1, 2, 3, 4, 5, 10, 15};
+    std::vector<f32> test_alphas = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, 30};
     std::transform(test_alphas.begin(), test_alphas.end(), test_alphas.begin(), to_radians);
 
     auto solvers = tiny::make_combination(meshes, backends);
@@ -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.02f || cd_rerr > 0.02f) return 1;
         }
     }
 

tests/test_non_linear.cpp

@@ -16,19 +16,23 @@ using namespace vlm;
 
 struct LiftCurveFunctor {
     virtual f32 operator()(f32 alpha) const = 0;
+    LiftCurveFunctor() = default;
+    virtual ~LiftCurveFunctor() = default;
 };
 
 struct SpallartLiftCurve : public LiftCurveFunctor {
     const f32 cl_0, a0, a1, cl_a0, cl_a1;
     SpallartLiftCurve(f32 cl_0_, f32 a0_, f32 a1_, f32 cl_a0_, f32 cl_a1_):
         cl_0(cl_0_), a0(a0_), a1(a1_), cl_a0(cl_a0_), cl_a1(cl_a1_) {}
+    ~SpallartLiftCurve() = default;
     inline f32 operator()(f32 alpha) const override {
         return cl_a0 * alpha + 0.5f * (cl_0 - cl_a1 * alpha) * (1.f + std::erf((alpha - a0) / a1));
     }
 };
 
 struct ThinAirfoilPolarLiftCurve : public LiftCurveFunctor{
     ThinAirfoilPolarLiftCurve() {}
+    ~ThinAirfoilPolarLiftCurve() = default;
     inline f32 operator()(f32 alpha) const override {
         return 2.f * PI_f * alpha;
     }

vlm/backends/cpu/src/vlm_backend_cpu.cpp

@@ -54,7 +54,7 @@ void BackendCPU::compute_lhs(const FlowData& flow) {
 
     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.panel_v.x.data(), m.panel_v.y.data(), m.panel_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()}
     };

vlm/backends/cpu/src/vlm_backend_cpu_kernels.ispc

@@ -87,7 +87,13 @@ export void kernel_influence(
     uniform float* uniform lhs,
     uniform uint64 ia, uniform uint64 lidx, uniform float sigma
     ) {
-    const uniform uint64 v0 = lidx + lidx / m.ns;
+    // const uniform uint64 v0 = lidx + lidx / m.ns;
+    // const uniform uint64 v1 = v0 + 1;
+    // const uniform uint64 v3 = v0 + m.ns+1;
+    // const uniform uint64 v2 = v3 + 1;
+    const uniform uint64 row = lidx / m.ns;
+    const uniform uint64 col = lidx % m.ns;
+    const uniform uint64 v0 = (2*row) * (m.ns+1) + col;
     const uniform uint64 v1 = v0 + 1;
     const uniform uint64 v3 = v0 + m.ns+1;
     const uniform uint64 v2 = v3 + 1;

vlm/dev/main.cpp

@@ -56,7 +56,7 @@ int main(int argc, char **argv) {
     });
 
     try {
-        // vlm::Executor::instance(1); // 1 thread
+        vlm::Executor::instance(1); // 1 thread
         LinearVLM solver(cfg);
         std::vector<f32> alphas = cfg().section("solver").get_vector<f32>("alphas", {0.0f});
         std::transform(alphas.begin(), alphas.end(), alphas.begin(), 

vlm/include/vlm_mesh.hpp

@@ -15,6 +15,7 @@ class Mesh {
     // ---------------------
     // All vertices stored in single SoA for result exporting
     // (stored in span major order)
+    SoA_3D_t<f32> panel_v; // size 2*(ncw)*(ns+1)
     SoA_3D_t<f32> v; // size (ncw+1)*(ns+1)
     // Offsets for indexing in connectivity array for each panel
     std::vector<u64> offsets = {}; // size nc*ns + 1

vlm/src/vlm_mesh.cpp

@@ -51,6 +51,7 @@ void Mesh::init() {
 void Mesh::alloc() {
     const u64 ncw = nc + nw;
     v.resize((ncw + 1) * (ns + 1));
+    panel_v.resize(2*ncw * (ns + 1));
     offsets.resize(nc * ns + 1);
     connectivity.resize(4 * nc * ns);
 
@@ -170,47 +171,64 @@ f32 Mesh::chord_length(const u64 j) const {
 }
 
 linalg::alias::float3 Mesh::get_v0(u64 i) const {
-    const u64 idx = i + i / ns;
-    return {v.x[idx], v.y[idx], v.z[idx]};
+    const u64 row = i / ns;
+    const u64 col = i % ns;
+    return {panel_v.x[(2*row+0) * (ns+1) + col], panel_v.y[(2*row+0) * (ns+1) + col], panel_v.z[(2*row+0) * (ns+1) + col]};
 }
 
 linalg::alias::float3 Mesh::get_v1(u64 i) const {
-    const u64 idx = i + i / ns + 1;
-    return {v.x[idx], v.y[idx], v.z[idx]};
+    const u64 row = i / ns;
+    const u64 col = i % ns;
+    return {panel_v.x[(2*row+0) * (ns+1) + col + 1], panel_v.y[(2*row+0) * (ns+1) + col + 1], panel_v.z[(2*row+0) * (ns+1) + col + 1]};
 }
 
 linalg::alias::float3 Mesh::get_v2(u64 i) const {
-    const u64 idx = i + i / ns + ns + 2;
-    return {v.x[idx], v.y[idx], v.z[idx]};
+    const u64 row = i / ns;
+    const u64 col = i % ns;
+    return {panel_v.x[(2*row+1) * (ns+1) + col + 1], panel_v.y[(2*row+1) * (ns+1) + col + 1], panel_v.z[(2*row+1) * (ns+1) + col + 1]};
 }
 
 linalg::alias::float3 Mesh::get_v3(u64 i) const {
-    const u64 idx = i + i / ns + ns + 1;
-    return {v.x[idx], v.y[idx], v.z[idx]};
+    const u64 row = i / ns;
+    const u64 col = i % ns;
+    return {panel_v.x[(2*row+1) * (ns+1) + col], panel_v.y[(2*row+1) * (ns+1) + col], panel_v.z[(2*row+1) * (ns+1) + col]};
 }
 
 void Mesh::update_wake(const linalg::alias::float3& freestream) {
-    const f32 chord_root = chord_length(0); 
-    const f32 off_x = freestream.x * 100.0f * chord_root;
-    const f32 off_y = freestream.y * 100.0f * chord_root;
-    const f32 off_z = freestream.z * 100.0f * chord_root;
+    const f32 chord_root = chord_length(0);
+    const linalg::alias::float3 vec_trans = linalg::normalize(freestream) * 100.0f * chord_root;
 
     const u64 v_ns = ns + 1;
     const u64 begin_trailing_edge = nb_vertices_wing()-v_ns;
     const u64 end_trailing_edge = nb_vertices_wing();
     // Add one layer of wake vertices
     // this can be parallelized (careful to false sharing tho)
     for (u64 i = begin_trailing_edge; i < end_trailing_edge; ++i) {
-        v.x[i + v_ns] = v.x[i] + off_x;
-        v.y[i + v_ns] = v.y[i] + off_y;
-        v.z[i + v_ns] = v.z[i] + off_z;
+        v.x[i + v_ns] = v.x[i] + vec_trans.x;
+        v.y[i + v_ns] = v.y[i] + vec_trans.y;
+        v.z[i + v_ns] = v.z[i] + vec_trans.z;
+    }
+
+    const u64 i = nc; // position on wake row
+    for (u64 j = 0; j < ns+1; j++) {
+        // copy trailing panels v2 and v3
+        panel_v.x[(2*i+0) * (ns+1) + j] = panel_v.x[(2*i-1) * (ns+1) + j];
+        panel_v.y[(2*i+0) * (ns+1) + j] = panel_v.y[(2*i-1) * (ns+1) + j];
+        panel_v.z[(2*i+0) * (ns+1) + j] = panel_v.z[(2*i-1) * (ns+1) + j];
+        // farfield wake vertices
+        panel_v.x[(2*i+1) * (ns+1) + j] = panel_v.x[(2*i+0) * (ns+1) + j] + vec_trans.x;
+        panel_v.y[(2*i+1) * (ns+1) + j] = panel_v.y[(2*i+0) * (ns+1) + j] + vec_trans.y;
+        panel_v.z[(2*i+1) * (ns+1) + j] = panel_v.z[(2*i+0) * (ns+1) + j] + vec_trans.z;
+
+        // std::printf("trailing: %f %f %f\n", panel_v.x[(2*i+0) * (ns+1) + j], panel_v.y[(2*i+0) * (ns+1) + j], panel_v.z[(2*i+0) * (ns+1) + j]);
+        // std::printf("wake: %f %f %f\n", panel_v.x[(2*i+1) * (ns+1) + j], panel_v.y[(2*i+1) * (ns+1) + j], panel_v.z[(2*i+1) * (ns+1) + j]);
     }
     compute_metrics_wake();
 }
 
 // https://publications.polymtl.ca/2555/1/2017_MatthieuParenteau.pdf (Eq 3.4 p21)
 void Mesh::correction_high_aoa(f32 alpha_rad) {
-    const f32 factor = 0.5f * alpha_rad / (std::sin(alpha_rad) + EPS); // correction factor
+    const f32 factor = (alpha_rad > EPS) ? 0.5f * alpha_rad / std::sin(alpha_rad) : 0.5f; // correction factor
     // Note: this can be vectorized and parallelized
     for (u64 i = 0; i < nb_panels_total(); i++) {
         // "chord vector" from center of leading line (v0-v1) to trailing line (v3-v2)
@@ -348,6 +366,34 @@ void Mesh::io_read_plot3d_structured(std::ifstream& f) {
             throw std::runtime_error("Mesh vertices should be ordered in chordwise direction");
         }
     }
+
+    // shift vertices according to quarter chord rule
+    for (u64 i = 0; i < nc; i++) {
+        for (u64 j = 0; j < ns+1; j++) {
+            linalg::alias::float3 pt0 = {v.x[i * nj + j], v.y[i * nj + j], v.z[i * nj + j]};
+            linalg::alias::float3 pt1 = {v.x[(i+1) * nj + j], v.y[(i+1) * nj + j], v.z[(i+1) * nj + j]};
+            linalg::alias::float3 vec_displacement = 0.25f * (pt1 - pt0);
+            pt0 += vec_displacement;
+            pt1 += vec_displacement;
+            panel_v.x[(2*i+0) * nj + j] = pt0.x;
+            panel_v.y[(2*i+0) * nj + j] = pt0.y;
+            panel_v.z[(2*i+0) * nj + j] = pt0.z;
+            panel_v.x[(2*i+1) * nj + j] = pt1.x;
+            panel_v.y[(2*i+1) * nj + j] = pt1.y;
+            panel_v.z[(2*i+1) * nj + j] = pt1.z;
+        }
+    }
+
+    // for (u64 i = 0; i < nc; i++) {
+    //     for (u64 j = 0; j < ns; j++) {
+    //         std::printf("panel: %llu\n", i*ns + j);
+    //         // in clockwise order starting from top left 
+    //         std::printf("v0: %f %f %f\n", panel_v.x[(2*i+0)*nj + j], panel_v.y[(2*i+0)*nj + j], panel_v.z[(2*i+0)*nj + j]);
+    //         std::printf("v1: %f %f %f\n", panel_v.x[(2*i+0)*nj + j + 1], panel_v.y[(2*i+0)*nj + j + 1], panel_v.z[(2*i+0)*nj + j + 1]);
+    //         std::printf("v2: %f %f %f\n", panel_v.x[(2*i+1)*nj + j + 1], panel_v.y[(2*i+1)*nj + j + 1], panel_v.z[(2*i+1)*nj + j + 1]);
+    //         std::printf("v3: %f %f %f\n", panel_v.x[(2*i+1)*nj + j], panel_v.y[(2*i+1)*nj + j], panel_v.z[(2*i+1)*nj + j]);
+    //     }
+    // }
 }
 
 void Mesh::io_read(const std::string& filename) {

xmake.lua

@@ -7,7 +7,7 @@ add_rules("mode.debug", "mode.release", "mode.releasedbg", "mode.asan")
 -- set_toolchains("cuda")
 
 -- set_toolset("cxx", "clang")
--- set_policy("build.sanitizer.address", true)
+set_policy("build.sanitizer.address", true)
 set_policy("build.warning", true)
 set_policy("build.cuda.devlink", true) -- magic
 set_policy("run.autobuild", true)