Skip to content

Commit

Permalink
more work
Browse files Browse the repository at this point in the history
  • Loading branch information
@samayala22
samayala22 committed Jun 19, 2024
1 parent 6e7db05 commit ca17c67
Show file tree
Hide file tree
Showing 12 changed files with 283 additions and 227 deletions.
12 changes: 8 additions & 4 deletions tests/test_vlm_elliptic_coeffs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ float compute_analytical_gamma(float y, float a, float b, float alpha) {
// return 0;
// }

int main(int argc, char **argv) {
int main(int /*argc*/, char ** /*argv*/) {
const float a = 1.0f; // wing chord root
const float b = 5.0f; // half wing span

Expand All @@ -115,11 +115,13 @@ int main(int argc, char **argv) {
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);

MeshGeom mesh_geom{};
auto solvers = tiny::make_combination(meshes, backends);
for (const auto& [mesh_name, backend_name] : solvers) {
mesh_io_read_file(mesh_name, &mesh_geom);
mesh_quarterchord(&mesh_geom);
LinearVLM solver{};
solver.mesh = create_mesh(mesh_name);
solver.backend = create_backend(backend_name, *solver.mesh);
solver.backend = create_backend(backend_name, &mesh_geom, 1);

for (u64 i = 0; i < test_alphas.size(); i++) {
const FlowData flow{test_alphas[i], 0.0f, 1.0f, 1.0f};
Expand All @@ -134,7 +136,9 @@ int main(int argc, char **argv) {
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.03f) return 1;
}
}

delete[] mesh_geom.vertices;
}
return 0;
}
6 changes: 5 additions & 1 deletion vlm/backends/cpu/include/vlm_backend_cpu.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@

namespace vlm {

class BackendCPU : public Backend {
class BackendCPU final : public Backend {
public:
BackendCPU(MeshGeom* mesh, u64 timesteps);
~BackendCPU();
Expand All @@ -24,6 +24,10 @@ class BackendCPU : public Backend {
void compute_delta_gamma() override;
void set_velocities(const linalg::alias::float3& vel) override;
void set_velocities(const SoA_3D_t<f32>& vels) override;

void mesh_metrics() override;
void mesh_update_wake(const linalg::alias::float3& freestream) override;
void mesh_correction_high_aoa(const f32 alpha) override;
private:
f32 mesh_mac(u64 j, u64 n) override; // mean chord
f32 mesh_area(const u64 i, const u64 j, const u64 m, const u64 n) override; // mean span
Expand Down
234 changes: 119 additions & 115 deletions vlm/backends/cpu/src/vlm_backend_cpu.cpp
View file

Large diffs are not rendered by default.

110 changes: 71 additions & 39 deletions vlm/backends/cpu/src/vlm_backend_cpu_kernels.ispc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,22 @@ inline uniform float length2(const uniform float3& v) {return LENGTH2(v);}
inline float length(const float3& v) {return LENGTH(v);}
inline uniform float length(const uniform float3& v) {return LENGTH(v);}

struct Mesh2 {
inline uniform float3 normalize(const uniform float3& v) {
uniform float l = length(v);
return {v.x/l, v.y/l, v.z/l};
}

inline float3 normalize(const float3& v) {
float l = length(v);
return {v.x/l, v.y/l, v.z/l};
}

template<typename F3>
inline F3 quad_normal(const F3& v0, const F3& v1, const F3& v2, const F3& v3) {
return normalize(cross(v3-v1, v2-v0));
}

export struct Mesh2 {
float* vertices; // (nc+nw+1)*(ns+1)*3
// TODO: evaluate if we really need those values or we can just compute on the fly in the kernels
float* normals; // nc*ns*3
Expand Down Expand Up @@ -145,6 +160,10 @@ export void kernel_wake_influence(

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

const uniform float* uniform vx = PTR_MESH_V(m, 0, 0, 0);
const uniform float* uniform vy = PTR_MESH_V(m, 0, 0, 1);
const uniform float* uniform vz = PTR_MESH_V(m, 0, 0, 2);

// Wake influence on wing
for (uint64 i = m->nc + m->nw - m->nwa; i < m->nc + m->nw; i++) {
foreach(j = 0 ... m->ns) {
Expand All @@ -153,10 +172,6 @@ export void kernel_wake_influence(
uint64 v2 = (i+1) * (m->ns+1) + j + 1;
uint64 v3 = (i+1) * (m->ns+1) + j;

const uniform float* uniform vx = PTR_MESH_V(m, 0, 0, 0);
const uniform float* uniform vy = PTR_MESH_V(m, 0, 0, 1);
const uniform float* uniform vz = PTR_MESH_V(m, 0, 0, 2);

// Loads
const float3 vertex0 = {vx[v0], vy[v0], vz[v0]};
const float3 vertex1 = {vx[v1], vy[v1], vz[v1]};
Expand Down Expand Up @@ -187,8 +202,11 @@ export void kernel_rollup(
uniform float* uniform rollup_vx, uniform float* uniform rollup_vy, uniform float* uniform rollup_vz,
uniform uint64 vidx, uniform float* uniform gamma, uniform float sigma
) {

const uniform float3 vertex_influenced = {m.v.x[vidx], m.v.y[vidx], m.v.z[vidx]};

const uniform float* uniform vx = PTR_MESH_V(m, 0, 0, 0);
const uniform float* uniform vy = PTR_MESH_V(m, 0, 0, 1);
const uniform float* uniform vz = PTR_MESH_V(m, 0, 0, 2);
const uniform float3 vertex_influenced = {vx[vidx], vy[vidx], vz[vidx]};

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

Expand All @@ -201,10 +219,10 @@ export void kernel_rollup(
uint64 v3 = (i+1) * (m->ns+1) + j;

// Loads
const float3 vertex0 = {m.v.x[v0], m.v.y[v0], m.v.z[v0]};
const float3 vertex1 = {m.v.x[v1], m.v.y[v1], m.v.z[v1]};
const float3 vertex2 = {m.v.x[v2], m.v.y[v2], m.v.z[v2]};
const float3 vertex3 = {m.v.x[v3], m.v.y[v3], m.v.z[v3]};
const float3 vertex0 = {vx[v0], vy[v0], vz[v0]};
const float3 vertex1 = {vx[v1], vy[v1], vz[v1]};
const float3 vertex2 = {vx[v2], vy[v2], vz[v2]};
const float3 vertex3 = {vx[v3], vy[v3], vz[v3]};

float3 ind = {0.0f, 0.0f, 0.0f};
kernel_symmetry(ind, vertex_influenced, vertex0, vertex1, sigma);
Expand All @@ -217,18 +235,18 @@ export void kernel_rollup(
}

// Wake influence
for (uint64 i = m->nc + m.nw - m.cw; i < m->nc + m.nw; i++) {
for (uint64 i = m->nc + m->nw - m->nwa; i < m->nc + m->nw; i++) {
foreach(j = 0 ... m->ns) {
uint64 v0 = (i+0) * (m->ns+1) + j;
uint64 v1 = (i+0) * (m->ns+1) + j + 1;
uint64 v2 = (i+1) * (m->ns+1) + j + 1;
uint64 v3 = (i+1) * (m->ns+1) + j;

// Loads
const float3 vertex0 = {m.v.x[v0], m.v.y[v0], m.v.z[v0]};
const float3 vertex1 = {m.v.x[v1], m.v.y[v1], m.v.z[v1]};
const float3 vertex2 = {m.v.x[v2], m.v.y[v2], m.v.z[v2]};
const float3 vertex3 = {m.v.x[v3], m.v.y[v3], m.v.z[v3]};
const float3 vertex0 = {vx[v0], vy[v0], vz[v0]};
const float3 vertex1 = {vx[v1], vy[v1], vz[v1]};
const float3 vertex2 = {vx[v2], vy[v2], vz[v2]};
const float3 vertex3 = {vx[v3], vy[v3], vz[v3]};

float3 ind = {0.0f, 0.0f, 0.0f};
kernel_symmetry(ind, vertex_influenced, vertex0, vertex1, sigma);
Expand All @@ -240,38 +258,49 @@ export void kernel_rollup(
}
}

rollup_vx[vidx] = m.v.x[vidx] + reduce_add(induced_vel.x) * dt;
rollup_vy[vidx] = m.v.y[vidx] + reduce_add(induced_vel.y) * dt;
rollup_vz[vidx] = m.v.z[vidx] + reduce_add(induced_vel.z) * dt;
rollup_vx[vidx] = vx[vidx] + reduce_add(induced_vel.x) * dt;
rollup_vy[vidx] = vy[vidx] + reduce_add(induced_vel.y) * dt;
rollup_vz[vidx] = vz[vidx] + reduce_add(induced_vel.z) * dt;
}

export uniform float kernel_trefftz_cd(
uniform const MeshProxy& m,
uniform const Mesh2* uniform m,
uniform float* uniform gamma,
uniform float* uniform trefftz_buffer,
uniform uint64 j, uniform uint64 n, uniform float sigma
) {
uniform uint64 begin = m.nb_panels + j;
const uniform uint64 nb_panels = m->ns * m->nc;
uniform uint64 begin = nb_panels + j;
uniform uint64 end = begin + n;
float cd = 0.0f;


const uniform float* uniform vx = PTR_MESH_V(m, 0, 0, 0);
const uniform float* uniform vy = PTR_MESH_V(m, 0, 0, 1);
const uniform float* uniform vz = PTR_MESH_V(m, 0, 0, 2);
const uniform float* uniform cx = PTR_MESH_C(m, 0, 0, 0);
const uniform float* uniform cy = PTR_MESH_C(m, 0, 0, 1);
const uniform float* uniform cz = PTR_MESH_C(m, 0, 0, 2);
const uniform float* uniform nx = PTR_MESH_N(m, 0, 0, 0);
const uniform float* uniform ny = PTR_MESH_N(m, 0, 0, 1);
const uniform float* uniform nz = PTR_MESH_N(m, 0, 0, 2);

// Compute the induced velocity of the streamwise wake vortex segments
for (uniform uint64 ia = m.nb_panels; ia < m.nb_panels + m.ns; ia++) {
const uniform uint64 v0 = ia + ia / m.ns;
for (uniform uint64 ia = nb_panels; ia < nb_panels + m->ns; ia++) {
const uniform uint64 v0 = ia + ia / m->ns;
const uniform uint64 v1 = v0 + 1;
const uniform uint64 v3 = v0 + m.ns+1;
const uniform uint64 v3 = v0 + m->ns+1;
const uniform uint64 v2 = v3 + 1;

// Broadcast vertices
const uniform float3 vertex0 = {m.v.x[v0], m.v.y[v0], m.v.z[v0]};
const uniform float3 vertex1 = {m.v.x[v1], m.v.y[v1], m.v.z[v1]};
const uniform float3 vertex2 = {m.v.x[v2], m.v.y[v2], m.v.z[v2]};
const uniform float3 vertex3 = {m.v.x[v3], m.v.y[v3], m.v.z[v3]};
const float3 vertex0 = {vx[v0], vy[v0], vz[v0]};
const float3 vertex1 = {vx[v1], vy[v1], vz[v1]};
const float3 vertex2 = {vx[v2], vy[v2], vz[v2]};
const float3 vertex3 = {vx[v3], vy[v3], vz[v3]};

uniform float gammaw = gamma[ia - m.ns];
uniform float gammaw = gamma[ia - m->ns];
foreach(ia2 = begin ... end) {
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]};
const float3 colloc = {cx[ia2], cy[ia2], cz[ia2]};
const float3 normal = {nx[ia2], ny[ia2], nz[ia2]};
float3 inf = {0.0f, 0.0f, 0.0f};

kernel_symmetry(inf, colloc, vertex1, vertex2, sigma);
Expand All @@ -281,25 +310,28 @@ export uniform float kernel_trefftz_cd(
}
// Perform the integration (Katz Plotkin, Low speed Aero | Eq 8.146)
foreach(i = j ... j + n) {
uint64 li = (m.nc-1) * m.ns + i;
float3 v0 = {m.v.x[i], m.v.y[i], m.v.z[i]};
float3 v1 = {m.v.x[i+1], m.v.y[i+1], m.v.z[i+1]};
uint64 li = (m->nc-1) * m->ns + i;
float3 v0 = {vx[i], vy[i], vz[i]};
float3 v1 = {vx[i+1], vy[i+1], vz[i+1]};
float dl = length(v1 - v0);
cd -= gamma[li] * trefftz_buffer[i] * dl; // used to have 0.5f * flow.rho
}
return reduce_add(cd);
}

export uniform float kernel_trefftz_cl(
uniform const MeshProxy& m,
uniform const Mesh2* uniform m,
uniform float* uniform gamma,
uniform uint64 j, uniform uint64 n
) {
float cl = 0.0f;
const uniform float* uniform vx = PTR_MESH_V(m, 0, 0, 0);
const uniform float* uniform vy = PTR_MESH_V(m, 0, 0, 1);
const uniform float* uniform vz = PTR_MESH_V(m, 0, 0, 2);
foreach(i = j ... j + n) {
uint64 li = (m.nc-1) * m.ns + i;
float3 v0 = {m.v.x[i], m.v.y[i], m.v.z[i]};
float3 v1 = {m.v.x[i+1], m.v.y[i+1], m.v.z[i+1]};
uint64 li = (m->nc-1) * m->ns + i;
float3 v0 = {vx[i], vy[i], vz[i]};
float3 v1 = {vx[i+1], vy[i+1], vz[i+1]};
float dl = length(v1 - v0);
cl += gamma[li]* dl; // used to have 0.5f * flow.rho
}
Expand Down
2 changes: 1 addition & 1 deletion vlm/backends/cuda/include/vlm_backend_cuda.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ struct MeshProxy {
SoA3D normal; // normals
};

class BackendCUDA : public Backend {
class BackendCUDA final : public Backend {
public:


Expand Down
1 change: 0 additions & 1 deletion vlm/include/vlm.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,6 @@ namespace vlm {

class Solver {
public:
std::unique_ptr<Mesh> mesh;
std::unique_ptr<Backend> backend;
Solver(const tiny::Config& cfg); // mesh filename & backend name
Solver() = default;
Expand Down
4 changes: 4 additions & 0 deletions vlm/include/vlm_backend.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,10 @@ class Backend {
virtual void set_velocities(const linalg::alias::float3& vel) = 0;
virtual void set_velocities(const SoA_3D_t<f32>& vels) = 0;

virtual void mesh_metrics() = 0;
virtual void mesh_update_wake(const linalg::alias::float3& freestream) = 0;
virtual void mesh_correction_high_aoa(const f32 alpha) = 0;

private:
virtual f32 mesh_mac(u64 j, u64 n) = 0; // mean chord
virtual f32 mesh_area(const u64 i, const u64 j, const u64 m, const u64 n) = 0; // mean span
Expand Down
3 changes: 2 additions & 1 deletion vlm/include/vlm_data.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,8 @@ namespace vlm {
struct Data {
f32* lhs = nullptr; // (ns*nc)^2
f32* rhs = nullptr; // ns*nc
f32* gamma = nullptr; // (nc+nw)*nw
f32* gamma = nullptr; // (nc+nw)*ns
f32* gamma_prev = nullptr; // nc*ns
f32* delta_gamma = nullptr; // nc*ns
f32* rollup_vertices = nullptr; // (nc+nw+1)*(ns+1)*3
f32* local_velocities = nullptr; // ns*nc*3
Expand Down
7 changes: 6 additions & 1 deletion vlm/include/vlm_mesh.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,10 @@
#include "tinyfwd.hpp"
#include "tinyview.hpp"

#define PTR_MESH_V(m, i,j,k) (m->vertices + j + i * (m->ns+1) + k * (m->nc+m->nw+1) * (m->ns+1))
#define PTR_MESH_C(m, i,j,k) (m->colloc + j + i * (m->ns) + k * (m->nc) * (m->ns))
#define PTR_MESH_N(m, i,j,k) (m->normals + j + i * (m->ns) + k * (m->nc) * (m->ns))

namespace vlm {

// C interface
Expand Down Expand Up @@ -41,7 +45,8 @@ struct MeshGeom {
u64 ns = 0;
};

void mesh_read_file(const std::string& filename, MeshGeom* mesh_geom);
void mesh_quarterchord(MeshGeom* mesh_geom);
void mesh_io_read_file(const std::string& filename, MeshGeom* mesh_geom);
void mesh_alloc(const malloc_f malloc, Mesh2* mesh, u64 nc, u64 ns, u64 nw, u64 nwa);
void mesh_free(const free_f free, Mesh2* mesh);
void mesh_copy(const memcpy_f memcpy, Mesh2* dst, const Mesh2* src);
Expand Down
Loading

0 comments on commit ca17c67

Please sign in to comment.