Merge pull request #6 from samayala22/ispc

Add ISPC backend

.github/workflows/linux.yaml

@@ -25,6 +25,12 @@ jobs:
       - name: Checkout repository
         uses: actions/checkout@v4
 
+      - name: Download ISPC
+        run: |
+          curl -L "https://github.com/ispc/ispc/releases/download/v1.22.0/ispc-v1.22.0-linux.tar.gz" -o ispc.tar.gz
+          tar -xzf ispc.tar.gz
+          echo "${PWD}/ispc-v1.22.0-linux/bin" >> $GITHUB_PATH
+
       # Install system dependencies (opengl)
       - name: Install system dependencies
         run: |

.github/workflows/windows.yaml

@@ -25,6 +25,12 @@ jobs:
       - name: Checkout repository
         uses: actions/checkout@v4
 
+      - name: Download ISPC
+        run: |
+          Invoke-WebRequest -Uri "https://github.com/ispc/ispc/releases/download/v1.22.0/ispc-v1.22.0-windows.zip" -OutFile "ispc.zip"
+          7z x ispc.zip
+          "ispc-v1.22.0-windows/bin" >> $env:GITHUB_PATH
+
       # Force xmake to a specific folder (for cache)
       - name: Set xmake env
         run: echo "XMAKE_GLOBALDIR=${{ runner.workspace }}/xmake-global" | Out-File -FilePath $env:GITHUB_ENV -Encoding utf8 -Append

config/elliptic.vlm

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

config/swept.vlm

@@ -7,7 +7,7 @@
 <solver>
 ref_point = [0.25, 0.0, 0.0]
 alphas = [5]
-backend = avx2
+backend = cpu
 
 <mesh>
 wake_included = false

headeronly/tinytimer.hpp

@@ -5,24 +5,24 @@
 
 namespace tiny {
 
-constexpr static long long us_in_s = 1'000'000LL;
-constexpr static long long us_in_ms = 1'000LL;
+constexpr static double us_in_s = 1'000'000;
+constexpr static double us_in_ms = 1'000;
 
-class Timer {
+class ScopedTimer {
 public:
-    Timer(const char* name) : m_name(name), m_start(std::chrono::high_resolution_clock::now()) {}
+    ScopedTimer(const char* name) : m_name(name), m_start(std::chrono::high_resolution_clock::now()) {}
 
-    ~Timer() {
+    ~ScopedTimer() {
         const auto m_end = std::chrono::high_resolution_clock::now();
-        const auto duration = static_cast<long long>(std::chrono::duration_cast<std::chrono::microseconds>(m_end - m_start).count());
+        const auto duration = static_cast<double>(std::chrono::duration_cast<std::chrono::microseconds>(m_end - m_start).count());
         std::printf("%s: ", m_name);
 
         if (duration > us_in_s) {
-            std::printf("%lld s\n", duration / us_in_s);
+            std::printf("%.2f s\n", duration / us_in_s);
         } else if (duration > us_in_ms) {
-            std::printf("%lld ms\n", duration / us_in_ms);
+            std::printf("%.0f ms\n", duration / us_in_ms);
         } else {
-            std::printf("%lld us\n", duration);
+            std::printf("%.0f us\n", duration);
         }
     }
 private:

mesh/infinite_rectangular_2x8.x

@@ -0,0 +1,23 @@
+       1
+       3        9        1
+     0.0000000000000000     0.50000000000000000      1.0000000000000000      0.0000000000000000
+    0.50000000000000000      1.0000000000000000      0.0000000000000000     0.50000000000000000
+     1.0000000000000000      0.0000000000000000     0.50000000000000000      1.0000000000000000
+     0.0000000000000000     0.50000000000000000      1.0000000000000000      0.0000000000000000
+    0.50000000000000000      1.0000000000000000      0.0000000000000000     0.50000000000000000
+     1.0000000000000000      0.0000000000000000     0.50000000000000000      1.0000000000000000
+     0.0000000000000000     0.50000000000000000      1.0000000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      3750.0000000000000
+     3750.0000000000005      3750.0000000000000      7500.0000000000000      7500.0000000000009
+     7500.0000000000000      11250.000000000000      11250.000000000000      11250.000000000000
+     15000.000000000000      15000.000000000002      15000.000000000000      18750.000000000000
+     18750.000000000000      18750.000000000000      17500.000000000000      17500.000000000000
+     17500.000000000000      13750.000000000000      13750.000000000000      13750.000000000000
+     10000.000000000000      10000.000000000000      10000.000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      0.0000000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      0.0000000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      0.0000000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      0.0000000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      0.0000000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      0.0000000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000

tests/test_linear.cpp

@@ -0,0 +1,141 @@
+#include <vector>
+#include <string>
+#include <fstream>
+#include <cmath>
+#include <cstdio>
+
+#include "tinycombination.hpp"
+
+#include "vlm.hpp"
+#include "vlm_utils.hpp"
+
+using namespace vlm;
+
+// area of whole wing
+float s_ref(float a, float b) {
+    return 0.5f * vlm::PI_f * a * b;
+}
+
+// wing aspect ratio
+float aspect_ratio(float a, float b) {
+    return (2*b)*(2*b) / s_ref(a, b);
+}
+
+float compute_analytical_cl(float alpha, float a, float b) {
+    return 2.0f * vlm::PI_f * alpha / (1.0f + 2.0f/aspect_ratio(a, b));
+}
+
+float compute_analytical_cd(float cl, float a, float b) {
+    return cl*cl / (vlm::PI_f * aspect_ratio(a, b));
+}
+
+float compute_analytical_gamma(float y, float a, float b, float alpha) {
+    const float gamma0 = compute_analytical_cl(alpha, a, b) * 1.0f * s_ref(a, b) / (vlm::PI_f * b);
+    const float ratio = y / b;
+    return gamma0 * std::sqrt(1.0f - ratio*ratio);
+}
+
+// bool elliptic_convergence() {
+//     tiny::Config cfg("../../../../config/elliptic.vlm");
+
+//     vlm::VLM vlm(cfg);
+
+//     std::vector<int> dimensions = {
+//         16, 32, 45, 64, 90, 128
+//     };
+
+//     const float a = 1.0f; // wing chord root
+//     const float b = 5.0f; // half wing span
+
+//     std::vector<double> norm_l1;
+//     std::vector<double> norm_l2;
+//     std::vector<double> norm_linf;
+
+//     const float alpha = 0.1f; // degrees
+
+//     for (const auto& dim : dimensions) {
+//         std::string filename = std::format("../../../../mesh/elliptic_{}x{}.x", dim, dim);
+//         vlm.mesh.io_read(filename);
+//         vlm.init();
+//         vlm::Solver solver(vlm.mesh, vlm.data, cfg);
+//         solver.run(alpha);
+
+//         double l1 = 0.0f;
+//         double l2 = 0.0f;
+//         double linf = 0.0f;
+//         int begin = (vlm.mesh.nc - 1) * vlm.mesh.ns;
+//         int end = vlm.mesh.nc * vlm.mesh.ns;
+//         // loop over last row of panels
+//         for (int i = begin; i < end; i++) {
+//             const float y = vlm.mesh.colloc.y[i];
+//             const float gamma = vlm.data.gamma[i];
+//             const float gamma_analytical = analytical_gamma(y, a, b, alpha);
+//             const double error = std::abs((gamma - gamma_analytical) / (gamma_analytical + 1e-7f));
+//             std::printf("y: %f, gamma: %f, gamma_analytical: %f, error: %f \n", y, gamma, gamma_analytical, error);
+
+//             l1 += error;
+//             l2 += error * error;
+//             linf = std::max(linf, error);
+//         }
+//         l1 /= (end - begin);
+//         l2 = std::sqrt(l2) / (end - begin);
+//         std::printf("L1: %f, L2: %f, Linf: %f\n", l1, l2, linf);
+
+//         norm_l1.push_back(l1);
+//         norm_l2.push_back(l2);
+//         norm_linf.push_back(linf);
+//     }
+
+//     double order_l1 = 0.0f;
+//     double order_l2 = 0.0f;
+//     double order_linf = 0.0f;
+
+//     auto order = [=](double norm0, double norm1, float dim0, float dim1) {
+//         return std::log(norm0 / norm1) / std::log((b/dim0)/(b/dim1));
+//     };
+
+//     for (int i = 0; i < dimensions.size() - 1; i++) {
+//         order_l1 += order(norm_l1[i], norm_l1[i+1], dimensions[i], dimensions[i+1]);
+//         order_l2 += order(norm_l2[i], norm_l2[i+1], dimensions[i], dimensions[i+1]);
+//         order_linf += order(norm_linf[i], norm_linf[i+1], dimensions[i], dimensions[i+1]);
+//     }
+//     order_l1 /= (dimensions.size() - 1);
+//     order_l2 /= (dimensions.size() - 1);
+//     order_linf /= (dimensions.size() - 1);
+//     std::printf("Order L1: %f, Order L2: %f, Order Linf: %f\n", order_l1, order_l2, order_linf);
+//     return 0;
+// }
+
+int main(int argc, char **argv) {
+    const float a = 1.0f; // wing chord root
+    const float b = 5.0f; // half wing span
+
+    const std::vector<std::string> meshes = {"../../../../mesh/elliptic_64x64.x"};
+    const std::vector<std::string> backends = {"cpu"};
+    std::vector<f32> test_alphas = {0, 1, 2, 3, 4, 5, 10, 15};
+    std::transform(test_alphas.begin(), test_alphas.end(), test_alphas.begin(), to_radians);
+
+    auto solvers = tiny::make_combination(meshes, backends);
+    for (const auto& [mesh_name, backend_name] : solvers) {
+        LinearVLM solver{};
+        solver.mesh = create_mesh(mesh_name);
+        solver.backend = create_backend(backend_name, *solver.mesh);
+
+        for (u64 i = 0; i < test_alphas.size(); i++) {
+            const FlowData flow{test_alphas[i], 0.0f, 1.0f, 1.0f};
+            const auto coeffs = solver.solve(flow);
+            const f32 analytical_cl = compute_analytical_cl(flow.alpha, a, b);
+            const f32 analytical_cd = compute_analytical_cd(analytical_cl, a, b);
+            const f32 cl_aerr = std::abs(coeffs.cl - analytical_cl);
+            const f32 cl_rerr = analytical_cl < EPS_f ? 0.f : cl_aerr / analytical_cl;
+            const f32 cd_aerr = std::abs(coeffs.cd - analytical_cd);
+            const f32 cd_rerr = analytical_cd < EPS_f ? 0.f : cd_aerr / analytical_cd;
+            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;
+        }
+    }
+
+    return 0;
+}

tests/test_non_linear.cpp

@@ -7,9 +7,11 @@
 #include <utility>
 
 #include "tinyinterpolate.hpp"
-#include "vlm.hpp"
 #include "tinycombination.hpp"
 
+#include "vlm.hpp"
+#include "vlm_utils.hpp"
+
 using namespace vlm;
 
 struct LiftCurveFunctor {
@@ -33,22 +35,14 @@ struct ThinAirfoilPolarLiftCurve : public LiftCurveFunctor{
 };
 
 template<typename T>
-void linspace(T start, T end, u32 n, std::vector<T>& out) {
+void linspace(T start, T end, u64 n, std::vector<T>& out) {
     out.resize(n);
     T step = (end - start) / (n - 1);
-    for (u32 i = 0; i < n; i++) {
+    for (u64 i = 0; i < n; i++) {
         out[i] = start + i * step;
     }
 }
 
-f32 to_radians(f32 degrees) {
-    return degrees * PI_f / 180.0f;
-}
-
-f32 to_degrees(f32 radians) {
-    return radians * 180.0f / PI_f;
-}
-
 template<typename T>
 void write_vector_pair(const std::string& filename, const std::vector<T>& vec1, const std::vector<T>& vec2) {
     assert(vec1.size() == vec2.size());
@@ -67,7 +61,7 @@ void write_vector_pair(const std::string& filename, const std::vector<T>& vec1,
 
 int main(int argc, char** argv) {
     const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_5x200.x"};
-    const std::vector<std::string> backends = {"avx2"};
+    const std::vector<std::string> backends = {"cpu"};
     std::vector<std::pair<std::string, std::unique_ptr<LiftCurveFunctor>>> lift_curves;
     lift_curves.emplace_back(std::make_pair("spallart1", std::make_unique<SpallartLiftCurve>(1.2f, 0.28f, 0.02f, 2.f*PI_f, 2.f*PI_f)));
     lift_curves.emplace_back(std::make_pair("spallart2", std::make_unique<SpallartLiftCurve>(0.72f, 0.28f, 0.04f, 2.f*PI_f, 1.5f*PI_f)));
@@ -106,16 +100,16 @@ int main(int argc, char** argv) {
             );
             db.profiles_pos.emplace_back(0.0f);
 
-            for (u32 i = 0; i < test_alphas.size(); i++) {
+            for (u64 i = 0; i < test_alphas.size(); i++) {
                 const FlowData flow{test_alphas[i], 0.0f, 1.0f, 1.0f};
                 auto coeffs = solver.solve(flow, db);
                 test_cl[i] = coeffs.cl;
                 std::printf(">>> Alpha: %.1f | CL = %.6f CD = %.6f CMx = %.6f CMy = %.6f CMz = %.6f\n", to_degrees(test_alphas[i]), coeffs.cl, coeffs.cd, coeffs.cm.x, coeffs.cm.y, coeffs.cm.z);
                 const f32 analytical_cl = (*lift_curve.second)(flow.alpha);
                 const f32 abs_error = std::abs(coeffs.cl - analytical_cl);
-                const f32 rel_error = 100.0f * abs_error / (analytical_cl + std::numeric_limits<f32>::epsilon());
-                std::printf(">>> Analytical: %.6f | Abs Error: %.3E | Relative Error: %.5f%% \n", analytical_cl, abs_error, rel_error);
-                if (rel_error > 1.f) return 1; // Failure
+                const f32 rel_error = abs_error / (analytical_cl + std::numeric_limits<f32>::epsilon());
+                std::printf(">>> Analytical: %.6f | Abs Error: %.3E | Relative Error: %.5f%% \n", analytical_cl, abs_error, rel_error*100.f);
+                if (rel_error > 0.01f) return 1; // Failure
             }
             write_vector_pair(lift_curve.first + "_nonlinear_cl", test_alphas, test_cl);
         }