diff --git a/data/theodorsen.py b/data/theodorsen.py index 1c7db41..8047351 100644 --- a/data/theodorsen.py +++ b/data/theodorsen.py @@ -25,8 +25,9 @@ def solve_ivp(x0: float, s0: float, sf: float, f: callable): # UVLM parameters rho = 1 # fluid density u_inf = 1 # freestream +ar = 4 # aspect ratio b = 0.5 # half chord -a = 10 # full span +a = ar / (2*b) # full span amplitudes = [0.1, 0.1, 0.1] reduced_frequencies = [0.5, 0.75, 1.5] @@ -76,9 +77,10 @@ def cl_theodorsen(t: float): uvlm_cl.append(float(cl)) axs["time"].plot(uvlm_t, uvlm_cl, label="UVLM (k=0.5)", linestyle="--") - axs["time"].plot(uvlm_t, uvlm_z, label="UVLM z (k=0.5)", linestyle="--") +axs["heave"].plot(uvlm_z[len(uvlm_cl)//2:], uvlm_cl[len(uvlm_cl)//2:], label="UVLM (k=0.5)", linestyle="--") + axs["time"].set_xlabel('t') axs["time"].set_ylabel('CL') axs["time"].grid(True, which='both', linestyle=':', linewidth=1.0, color='gray') diff --git a/mesh/rectangular_5x10.x b/mesh/rectangular_5x10.x new file mode 100644 index 0000000..9102ef5 --- /dev/null +++ b/mesh/rectangular_5x10.x @@ -0,0 +1,53 @@ + 1 + 6 11 1 + 0.0000000000000000 0.19999999999999996 0.39999999999999991 0.59999999999999998 + 0.80000000000000004 1.0000000000000000 0.0000000000000000 0.19999999999999998 + 0.39999999999999991 0.59999999999999998 0.80000000000000016 1.0000000000000000 + 0.0000000000000000 0.20000000000000001 0.39999999999999991 0.59999999999999998 + 0.79999999999999993 1.0000000000000000 0.0000000000000000 0.20000000000000004 + 0.39999999999999991 0.59999999999999998 0.79999999999999993 1.0000000000000000 + 0.0000000000000000 0.20000000000000001 0.39999999999999991 0.59999999999999998 + 0.80000000000000027 1.0000000000000000 0.0000000000000000 0.20000000000000007 + 0.39999999999999991 0.60000000000000009 0.79999999999999993 1.0000000000000000 + 0.0000000000000000 0.20000000000000009 0.39999999999999991 0.60000000000000009 + 0.79999999999999993 1.0000000000000000 0.0000000000000000 0.20000000000000012 + 0.39999999999999991 0.60000000000000009 0.80000000000000004 1.0000000000000000 + 0.0000000000000000 0.20000000000000012 0.39999999999999991 0.60000000000000009 + 0.79999999999999982 1.0000000000000000 0.0000000000000000 0.20000000000000015 + 0.39999999999999991 0.60000000000000009 0.79999999999999982 1.0000000000000000 + 0.0000000000000000 0.20000000000000018 0.39999999999999991 0.60000000000000009 + 0.79999999999999982 1.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 0.40000000000000002 0.39999999999999986 + 0.39999999999999947 0.39999999999999913 0.39999999999999897 0.39999999999999858 + 0.80000000000000004 0.79999999999999971 0.79999999999999893 0.79999999999999827 + 0.79999999999999793 0.79999999999999716 1.2000000000000002 1.2000000000000004 + 1.1999999999999997 1.1999999999999997 1.1999999999999991 1.1999999999999993 + 1.6000000000000001 1.6000000000000003 1.6000000000000005 1.6000000000000010 + 1.6000000000000012 1.6000000000000014 2.0000000000000000 2.0000000000000000 + 2.0000000000000000 2.0000000000000000 2.0000000000000000 2.0000000000000000 + 2.4000000000000004 2.3999999999999999 2.3999999999999995 2.3999999999999995 + 2.3999999999999981 2.3999999999999986 2.8000000000000003 2.8000000000000003 + 2.8000000000000012 2.8000000000000007 2.8000000000000007 2.8000000000000007 + 3.2000000000000002 3.2000000000000002 3.1999999999999997 3.1999999999999997 + 3.1999999999999993 3.1999999999999993 3.5999999999999996 3.5999999999999992 + 3.6000000000000005 3.6000000000000005 3.6000000000000010 3.6000000000000014 + 4.0000000000000000 4.0000000000000000 4.0000000000000000 4.0000000000000000 + 4.0000000000000000 4.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 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 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 diff --git a/tests/test_uvlm_theodorsen.cpp b/tests/test_uvlm_theodorsen.cpp index 404e068..16f2d02 100644 --- a/tests/test_uvlm_theodorsen.cpp +++ b/tests/test_uvlm_theodorsen.cpp @@ -64,13 +64,12 @@ void dump_buffer(std::ofstream& stream, T* start, T* end) { int main() { // const std::vector meshes = {"../../../../mesh/infinite_rectangular_2x8.x"}; - const std::vector meshes = {"../../../../mesh/infinite_rectangular_5x200.x"}; + const std::vector meshes = {"../../../../mesh/rectangular_5x10.x"}; const std::vector backends = get_available_backends(); auto solvers = tiny::make_combination(meshes, backends); - // Geometry const f32 b = 0.5f; // half chord @@ -78,7 +77,7 @@ int main() { const f32 t_final = 30.0f; const f32 u_inf = 1.0f; // freestream velocity const f32 amplitude = 0.1f; // amplitude of the wing motion - const f32 k = 0.5f; // reduced frequency + const f32 k = 1.5f; // reduced frequency const f32 omega = k * u_inf / (2*b); @@ -165,6 +164,8 @@ int main() { std::cout << "position: " << base_vertex << "\n"; const FlowData flow{linalg::alias::float3{-base_velocity[0], -base_velocity[1], -base_velocity[2]}, 1.0f}; + std::cout << "lift axis: " << flow.lift_axis << "\n"; + backend->compute_rhs(flow); backend->add_wake_influence(flow); backend->lu_solve(); diff --git a/vlm/backends/cpu/src/vlm_backend_cpu.cpp b/vlm/backends/cpu/src/vlm_backend_cpu.cpp index 4009c56..37e8a7f 100644 --- a/vlm/backends/cpu/src/vlm_backend_cpu.cpp +++ b/vlm/backends/cpu/src/vlm_backend_cpu.cpp @@ -90,7 +90,7 @@ void BackendCPU::compute_lhs(const FlowData& flow) { init.precede(wing_pass, cond); wing_pass.precede(sync); - cond.precede(wake_pass, sync); + cond.precede(wake_pass, sync); // 0 and 1 wake_pass.precede(back); back.precede(cond); @@ -125,18 +125,49 @@ void BackendCPU::add_wake_influence(const FlowData& flow) { }; // loop over wake rows - for (u64 i = 0; i < mesh.current_nw; i++) { - const u64 wake_row_start = (m.nc + m.nw - i - 1) * m.ns; - std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out - // Actually fill the wake buffer - // parallel for - for (u64 j = 0; j < mesh.ns; j++) { // loop over columns - const u64 lidx = wake_row_start + j; // TODO: replace this ASAP - ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas); - } + // for (u64 i = 0; i < mesh.current_nw; i++) { + // const u64 wake_row_start = (m.nc + m.nw - i - 1) * m.ns; + // std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out + // // Actually fill the wake buffer + // // parallel for + // for (u64 j = 0; j < mesh.ns; j++) { // loop over columns + // const u64 lidx = wake_row_start + j; // TODO: replace this ASAP + // ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas); + // } + + // cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, -1.0f, wake_buffer.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, rhs.data(), 1); + // } + + u64 idx = 0; + u64 wake_row_start = (m.nc + m.nw - 1) * m.ns; + auto init = taskflow.placeholder(); + auto cond = taskflow.emplace([&]{ + return idx < mesh.current_nw ? 0 : 1; // 0 means continue, 1 means break + }); + auto zero_buffer = taskflow.emplace([&]{ + std::fill(wake_buffer.begin(), wake_buffer.end(), 0.0f); // zero out + }); + auto wake_influence = taskflow.for_each_index((u64)0, m.ns, [&] (u64 j) { + const u64 lidx = wake_row_start + j; + ispc::kernel_influence(mesh_proxy, wake_buffer.data(), j, lidx, flow.sigma_vatistas); + }); + auto back = taskflow.emplace([&]{ cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, -1.0f, wake_buffer.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, rhs.data(), 1); - } + + idx++; + wake_row_start -= m.ns; + return 0; // 0 means continue + }); + auto sync = taskflow.placeholder(); + + init.precede(cond); + cond.precede(zero_buffer, sync); + zero_buffer.precede(wake_influence); + wake_influence.precede(back); + back.precede(cond); + + Executor::get().run(taskflow).wait(); } void BackendCPU::shed_gamma() { @@ -216,9 +247,16 @@ f32 BackendCPU::compute_coefficient_unsteady_cl(const FlowData& flow, f32 dt, co // Steady part: const linalg::alias::float3 v0 = mesh.get_v0(li); const linalg::alias::float3 v1 = mesh.get_v1(li); + const linalg::alias::float3 v2 = mesh.get_v2(li); + const linalg::alias::float3 v3 = mesh.get_v3(li); + + linalg::alias::float3 force_steady = {0.0f, 0.0f, 0.0f}; + force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v1 - v0); + force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v2 - v1); + force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v3 - v2); + force_steady += flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, v0 - v3); + // Leading edge vector pointing outward from wing root - const linalg::alias::float3 dl = v1 - v0; - const linalg::alias::float3 force_steady = flow.rho * delta_gamma[li] * linalg::cross(flow.freestream, dl); cl += linalg::dot(force_steady, flow.lift_axis); // Unsteady part (Simpson method)