hardstuck

data/debug.py

@@ -0,0 +1,11 @@
+import numpy as np
+
+lhs = np.array([-1.27324,3.81972e-10,0.424413,3.81972e-10,3.81972e-10,-1.27324,-1.46841e-08,0.424413,0.424413,3.81972e-10,-1.27324,3.81972e-10,1.54481e-08,0.424413,3.81972e-10,-1.27324]).reshape((4,4)).T
+
+rhs = np.array([-0.0871557519] * 4)
+wake_buf = np.array([0.0848438,3.80518e-10,0.424593,3.80518e-10,6.40765e-09,0.0848438,1.54434e-08,0.424593]).reshape((2,4)).T
+gamma = np.array([0.102677956, 0.102677956])
+
+
+print(wake_buf)
+print(rhs - wake_buf @ gamma)

data/theodorsen.py

@@ -14,9 +14,6 @@ def derivative2(f, x):
 def solve_ivp(x0: float, s0: float, sf: float, f: callable):
     return spi.solve_ivp(f, [s0, sf], [x0]).y[-1] # return only the result at t=sf
 
-def pade_approximation(k: float):
-    return (0.5177*k*k + 0.2752*k + 0.01576) / (k*k + 0.3414*k + 0.01582)
-
 # Some info in Katz Plotkin p414 (eq 13.73a)
 # Jone's approximation of Wagner function
 b0 = 1
@@ -51,8 +48,8 @@ def atime(t: float): return 2. * u_inf * t / c
 
 for amp, k in zip(amplitudes, reduced_frequencies):
     # heave parameters
-    amplitude = amp / (2*b)
-    omega = k * u_inf / (2*b) # pitch frequency
+    amplitude = amp / c
+    omega = k * u_inf / c # pitch frequency
 
     # sudden acceleration
     # def pitch(t): return np.radians(5)
@@ -78,14 +75,9 @@ def cl_theodorsen(t: float): # using Wagner functions and Kholodar formulation
         L_m = rho * b * b * np.pi * (u_inf * derivative(pitch, t) + derivative2(heave, t) - b * pitch_axis * derivative2(pitch, t))
         L_c = -2 * np.pi * rho * u_inf * b * (-(b0 + b1 + b2) * w(t) + x1(t) + x2(t))
         return (L_m + L_c) / (0.5 * rho * u_inf * u_inf * c)
-
-    # def cl_theodorsen(t: float): # using Pade approximation
-        # return 0.5 * np.pi * (derivative2(heave, t) + derivative(pitch, t) - 0.5 * pitch_axis * derivative2(pitch, t)) + 2.0 * np.pi * (pitch(t) + derivative(heave, t) + 0.5 * derivative(pitch, t) * (0.5 - pitch_axis)) * pade_approximation(k)
-        # L = rho * b * b * np.pi * (u_inf * derivative(pitch, t) + derivative2(heave, t) - b * pitch_axis * derivative2(pitch, t)) + 2 * np.pi * rho * u_inf * b * (u_inf * pitch(t) + derivative(heave, t) + b * (0.5 - pitch_axis) * derivative(pitch, t)) * pade_approximation(k)
-        # return L / (0.5 * rho * u_inf * u_inf * a * (2*b))
 
     cl_theo = np.array([cl_theodorsen(ti) for ti in vec_t])
-    coord_z = np.array([-heave(ti) / (2*b) for ti in vec_t])
+    coord_z = np.array([-heave(ti) / c for ti in vec_t])
 
     axs["time"].plot(vec_t, cl_theo, label=f"k={k}")
     axs["heave"].plot(coord_z[int(nb_pts//2):], cl_theo[int(nb_pts//2):], label=f"k={k}")
@@ -104,16 +96,13 @@ def cl_theodorsen(t: float): # using Wagner functions and Kholodar formulation
 uvlm_cl = np.array(uvlm_cl)
 analytical_cl = np.array([np.interp(ut, vec_t, cl_theo) for ut in uvlm_t])
 err_rel = np.abs((uvlm_cl-analytical_cl)/analytical_cl)
-quotient = np.abs(uvlm_cl / analytical_cl)
 print("Avg rel error: ", 100.0 * np.mean(err_rel))
-print("Quotient: ", np.mean(quotient))
 
-axs["error"].plot(uvlm_t, err_rel, label="rel")
-axs["error"].plot(uvlm_t, quotient, label="quotient")
+axs["error"].plot(uvlm_t, 100.0 * err_rel, label="rel (%)")
 axs["time"].scatter(uvlm_t, uvlm_cl, label="UVLM (k=0.5)", facecolors='none', edgecolors='b', s=15)
 axs["heave"].scatter(uvlm_z[len(uvlm_cl)//4:], uvlm_cl[len(uvlm_cl)//4:], label="UVLM (k=0.5)", facecolors='none', edgecolors='b', s=15)
 
-# axs["time"].plot(t, [0.548311] * len(t), label="VLM (alpha=5)")
+# axs["time"].plot(vec_t, [0.548311] * len(vec_t), label="VLM (alpha=5)")
 
 axs["time"].set_xlabel('t')
 axs["time"].set_ylabel('CL')
@@ -131,4 +120,6 @@ def cl_theodorsen(t: float): # using Wagner functions and Kholodar formulation
 axs["heave"].legend()
 
 plt.suptitle("Verification of UVLM with Theodorsen harmonic heave motion")
+# plt.suptitle("Verification of UVLM with sudden acceleration motion")
+
 plt.show()

mesh/infinite_rectangular_2x2.x

@@ -0,0 +1,11 @@
+       1
+       3        3        1
+     0.0000000000000000     0.50000000000000000      1.0000000000000000      0.0000000000000000
+    0.50000000000000000      1.0000000000000000      0.0000000000000000     0.50000000000000000
+     1.0000000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      20000.000000000000
+     20000.000000000000      20000.000000000000      40000.000000000000      40000.000000000000
+     40000.000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      0.0000000000000000
+     0.0000000000000000      0.0000000000000000      0.0000000000000000      0.0000000000000000
+     0.0000000000000000

tests/test_uvlm_theodorsen.cpp

@@ -62,9 +62,18 @@ void dump_buffer(std::ofstream& stream, T* start, T* end) {
     stream << "\n";
 }
 
+template<typename T>
+void print_buffer(const T* start, u64 size) {
+    std::cout << "[";
+    for (u64 i = 0; i < size; i++) {
+        std::cout << start[i] << ",";
+    }
+    std::cout << "]\n";
+}
+
 int main() {
     // const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_2x8.x"};
-    const std::vector<std::string> meshes = {"../../../../mesh/long_rectangular_10x10.x"};
+    const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_5x20.x"};
 
     const std::vector<std::string> backends = get_available_backends();
 
@@ -120,7 +129,7 @@ int main() {
         std::ofstream wing_data("wing_data.txt");
         std::ofstream wake_data("wake_data.txt");
         std::ofstream cl_data("cl_data.txt");
-        
+
         wing_data << mesh->nc << " " << mesh->ns << "\n";
         wing_data << vec_t.size() - 1 << "\n\n";
 
@@ -171,7 +180,7 @@ int main() {
             if (i > 0) {
                 // TODO: this should take a vector of local velocities magnitude because it can be different for each point on the mesh
                 const f32 cl_unsteady = backend->compute_coefficient_unsteady_cl(flow, dt, mesh->s_ref, 0, mesh->ns);
-                // const f32 cl_unsteady = backend->compute_coefficient_cl(flow);
+
                 std::printf("t: %f, CL: %f\n", t, cl_unsteady);
                 #ifdef DEBUG_DISPLACEMENT_DATA
                 cl_data << t << " " << mesh->v.z[0] << " " << cl_unsteady << "\n";
@@ -181,6 +190,6 @@ int main() {
             mesh->move(kinematics.relative_displacement(t, t+dt));
         }
     }
-  
+
     return 0;
 }

vlm/backends/cpu/src/vlm_backend_cpu.cpp

@@ -19,6 +19,16 @@
 #include <cblas.h>
 
 using namespace vlm;
+
+template<typename T>
+void print_buffer(const T* start, u64 size) {
+    std::cout << "[";
+    for (u64 i = 0; i < size; i++) {
+        std::cout << start[i] << ",";
+    }
+    std::cout << "]\n";
+}
+
 using namespace linalg::ostream_overloads;
 
 BackendCPU::~BackendCPU() = default; // Destructor definition