Skip to content

Commit

Permalink
compute induced wake velocities
Browse files Browse the repository at this point in the history 
disgustingly inefficient implementation
  • Loading branch information
@samayala22
samayala22 committed May 1, 2024
1 parent b9ffab3 commit 1422de4
Show file tree
Hide file tree
Showing 5 changed files with 139 additions and 56 deletions.
13 changes: 7 additions & 6 deletions data/theodorsen.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ def atime(t: float): return 2. * u_inf * t / c

cycles = 4
amplitudes = [0.1]
reduced_frequencies = [0.75]
reduced_frequencies = [0.5]

t_final = cycles * 2 * np.pi / (reduced_frequencies[0] * 2.0 * u_inf / c) # 4 cycles
nb_pts = 500
Expand All @@ -57,20 +57,20 @@ def atime(t: float): return 2. * u_inf * t / c

for amp, k in zip(amplitudes, reduced_frequencies):
# heave parameters
amplitude = amp / c
amplitude = amp
omega = k * 2.0 * u_inf / c # pitch frequency

# sudden acceleration
# def pitch(t): return np.radians(5)
# 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)
Expand Down Expand Up @@ -119,6 +119,7 @@ def cl_theodorsen(t: float): # using Wagner functions and Kholodar formulation
print(f"Error: {100 * error / np.max(np.abs(analytical_cl)):.3f}%", )

axs["time"].scatter(uvlm_t, uvlm_cl, label="UVLM", facecolors='none', edgecolors='r', s=20)
# axs["time"].scatter(uvlm_t, uvlm_z, label="UVLM z", facecolors='none', edgecolors='r', s=20)
axs["heave"].scatter(uvlm_z[uvlm_cycle_idx:], uvlm_cl[uvlm_cycle_idx:], label="UVLM", facecolors='none', edgecolors='r', s=20)

# axs["time"].plot(vec_t, [0.548311] * len(vec_t), label="VLM (alpha=5)")
Expand Down
7 changes: 4 additions & 3 deletions tests/test_uvlm_theodorsen.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_10x5.x"};
const std::vector<std::string> meshes = {"../../../../mesh/infinite_rectangular_2x2.x"};

const std::vector<std::string> backends = get_available_backends();

Expand All @@ -90,9 +90,10 @@ int main() {
const f32 cycles = 4.0f;
const f32 u_inf = 1.0f; // freestream velocity
const f32 amplitude = 0.1f; // amplitude of the wing motion
const f32 k = 0.75f; // reduced frequency
const f32 k = 0.5; // reduced frequency
const f32 omega = k * 2.0f * u_inf / (2*b);
const f32 t_final = cycles * 2.0f * PI_f / omega; // 4 periods
//const f32 t_final = 5.0f;

Kinematics kinematics{};

Expand Down Expand Up @@ -153,7 +154,7 @@ int main() {
f32 dt = segment_chord / kinematics.velocity_magnitude(t, {trailing_vertices.x[0], trailing_vertices.y[0], trailing_vertices.z[0], 1.0f});
for (u64 i = 1; i < trailing_vertices.size; i++) {
dt = std::min(dt, segment_chord / kinematics.velocity_magnitude(t, {trailing_vertices.x[i], trailing_vertices.y[i], trailing_vertices.z[i], 1.0f}));
}
}

auto transform = kinematics.relative_displacement(t, t+dt);
for (u64 i = 0; i < trailing_vertices.size; i++) {
Expand Down
6 changes: 6 additions & 0 deletions vlm/backends/cpu/include/vlm_backend_cpu.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,12 @@ namespace vlm {
class BackendCPU : public Backend {
public:
std::vector<f32> lhs;
std::vector<f32> uw;
std::vector<f32> vw;
std::vector<f32> ww;
std::vector<f32> panel_uw;
std::vector<f32> panel_vw;
std::vector<f32> panel_ww;
std::vector<f32> wake_buffer; // (ns*nc) * ns
std::vector<f32> rhs;
std::vector<i32> ipiv;
Expand Down
124 changes: 77 additions & 47 deletions vlm/backends/cpu/src/vlm_backend_cpu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,15 @@ BackendCPU::~BackendCPU() = default; // Destructor definition
BackendCPU::BackendCPU(Mesh& mesh) : Backend(mesh) {
lhs.resize(mesh.nb_panels_wing() * mesh.nb_panels_wing());
wake_buffer.resize(mesh.nb_panels_wing() * mesh.ns);

uw.resize(mesh.nb_panels_wing() * mesh.ns);
vw.resize(mesh.nb_panels_wing() * mesh.ns);
ww.resize(mesh.nb_panels_wing() * mesh.ns);

panel_uw.resize(mesh.nb_panels_wing());
panel_vw.resize(mesh.nb_panels_wing());
panel_ww.resize(mesh.nb_panels_wing());

rhs.resize(mesh.nb_panels_wing());
ipiv.resize(mesh.nb_panels_wing());
gamma.resize((mesh.nc + mesh.nw) * mesh.ns); // store wake gamma as well
Expand Down Expand Up @@ -136,50 +145,60 @@ void BackendCPU::add_wake_influence() {
{m.normal.x.data(), m.normal.y.data(), m.normal.z.data()}
};

std::fill(panel_uw.begin(), panel_uw.end(), 0.0f);
std::fill(panel_vw.begin(), panel_vw.end(), 0.0f);
std::fill(panel_ww.begin(), panel_ww.end(), 0.0f);

// 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, 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
std::fill(uw.begin(), uw.end(), 0.0f);
std::fill(vw.begin(), vw.end(), 0.0f);
std::fill(ww.begin(), ww.end(), 0.0f);
// 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_influence2(mesh_proxy, wake_buffer.data(), uw.data(), vw.data(), ww.data(), j, lidx, 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);
// }
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);
cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, 1.0f, uw.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, panel_uw.data(), 1);
cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, 1.0f, vw.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, panel_vw.data(), 1);
cblas_sgemv(CblasColMajor, CblasNoTrans, m.nb_panels_wing(), m.ns, 1.0f, ww.data(), m.nb_panels_wing(), gamma.data() + wake_row_start, 1, 1.0f, panel_ww.data(), 1);
}

u64 idx = 0;
u64 wake_row_start = (m.nc + m.nw - 1) * m.ns;
// 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, 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);
// 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, 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();
// 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);
// 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();
// Executor::get().run(taskflow).wait();
}

void BackendCPU::shed_gamma() {
Expand Down Expand Up @@ -249,11 +268,22 @@ f32 BackendCPU::compute_coefficient_unsteady_cl(const SoA_3D_t<f32>& vel, f32 dt
f32 cl = 0.0f;
const f32 rho = 1.0f; // TODO: remove hardcoded rho

// for (auto& elem : rhs) {
// std::cout << elem << " ";
// }
// std::cout << "\n";

for (u64 u = 0; u < mesh.nc; u++) {
for (u64 v = j; v < j + n; v++) {
const u64 li = u * mesh.ns + v; // linear index

const linalg::alias::float3 freestream{vel.x[li], vel.y[li], vel.z[li]};
linalg::alias::float3 freestream{vel.x[li], vel.y[li], vel.z[li]};
//std::cout << "Without: " << freestream << "\n";
freestream.x += panel_uw[li];
freestream.y += panel_vw[li];
freestream.z += panel_ww[li];
// std::cout << "With: " << freestream << "\n";

const linalg::alias::float3 lift_axis = compute_lift_axis(freestream);

const linalg::alias::float3 v0 = mesh.get_v0(li);
Expand All @@ -266,17 +296,17 @@ f32 BackendCPU::compute_coefficient_unsteady_cl(const SoA_3D_t<f32>& vel, f32 dt
const f32 gamma_dt = (gamma[li] - gamma_prev[li]) / dt; // backward difference

// Joukowski method
force += rho * delta_gamma[li] * linalg::cross(freestream, v1 - v0);
force += rho * gamma_dt * mesh.area[li] * normal;
// force += rho * delta_gamma[li] * linalg::cross(freestream, v1 - v0);
// force += rho * gamma_dt * mesh.area[li] * normal;

// Katz Plotkin method
// linalg::alias::float3 delta_p = {0.0f, 0.0f, 0.0f};
// const f32 delta_gamma_i = (u == 0) ? gamma[li] : gamma[li] - gamma[(u-1) * mesh.ns + v];
// const f32 delta_gamma_j = (v == 0) ? gamma[li] : gamma[li] - gamma[u * mesh.ns + v - 1];
// delta_p += rho * linalg::dot(freestream, linalg::normalize(v1 - v0)) * delta_gamma_j / mesh.panel_width_y(u, v);
// delta_p += rho * linalg::dot(freestream, linalg::normalize(v3 - v0)) * delta_gamma_i / mesh.panel_length(u, v);
// delta_p += gamma_dt / dt;
// force = (delta_p * mesh.area[li]) * normal;
linalg::alias::float3 delta_p = {0.0f, 0.0f, 0.0f};
const f32 delta_gamma_i = (u == 0) ? gamma[li] : gamma[li] - gamma[(u-1) * mesh.ns + v];
const f32 delta_gamma_j = (v == 0) ? gamma[li] : gamma[li] - gamma[u * mesh.ns + v - 1];
delta_p += rho * linalg::dot(freestream, linalg::normalize(v1 - v0)) * delta_gamma_j / mesh.panel_width_y(u, v);
delta_p += rho * linalg::dot(freestream, linalg::normalize(v3 - v0)) * delta_gamma_i / mesh.panel_length(u, v);
delta_p += gamma_dt;
force = (delta_p * mesh.area[li]) * normal;

// force /= linalg::length2(freestream);

Expand Down
45 changes: 45 additions & 0 deletions vlm/backends/cpu/src/vlm_backend_cpu_kernels.ispc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,51 @@ export void kernel_influence(
}
}

export void kernel_influence2(
uniform const MeshProxy& m,
uniform float* uniform lhs, uniform float* uniform uw, uniform float* uniform vw, uniform float* uniform ww,
uniform uint64 ia, uniform uint64 lidx, uniform float sigma
) {

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;


// 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]};
uniform float3 vertex2 = {m.v.x[v2], m.v.y[v2], m.v.z[v2]};
uniform float3 vertex3 = {m.v.x[v3], m.v.y[v3], m.v.z[v3]};

// print("Influencer %: \n", lidx);
// print("Vertex 0: % % %\n", vertex0.x, vertex0.y, vertex0.z);
// print("Vertex 1: % % %\n", vertex1.x, vertex1.y, vertex1.z);
// print("Vertex 2: % % %\n", vertex2.x, vertex2.y, vertex2.z);
// print("Vertex 3: % % %\n", vertex3.x, vertex3.y, vertex3.z);

foreach(ia2 = 0 ... m.nb_panels) {
const float3 colloc = {m.colloc.x[ia2], m.colloc.y[ia2], m.colloc.z[ia2]};
const float3 normal = {m.normal.x[ia2], m.normal.y[ia2], m.normal.z[ia2]};

// print("Influenced: %\n", ia2);
// print("Colloc: \n % \n % \n %\n", colloc.x, colloc.y, colloc.z);
// print("Normal: \n % \n % \n %\n", normal.x, normal.y, normal.z);

float3 inf = {0.0f, 0.0f, 0.0f};

kernel_symmetry(inf, colloc, vertex0, vertex1, sigma);
kernel_symmetry(inf, colloc, vertex1, vertex2, sigma);
kernel_symmetry(inf, colloc, vertex2, vertex3, sigma);
kernel_symmetry(inf, colloc, vertex3, vertex0, sigma);
lhs[ia * m.nb_panels + ia2] += dot(inf, normal); // store
uw[ia * m.nb_panels + ia2] += inf.x;
vw[ia * m.nb_panels + ia2] += inf.y;
ww[ia * m.nb_panels + ia2] += inf.z;
}
}

export uniform float kernel_trefftz_cd(
uniform const MeshProxy& m,
uniform float* uniform gamma,
Expand Down

0 comments on commit 1422de4

Please sign in to comment.