initial commit

vlm/backends/cpu/src/vlm_backend_cpu_kernels.ispc

@@ -35,9 +35,9 @@ export struct SoA3D {
 };
 
 export struct MeshProxy {
-    uniform unsigned int64 ns;
-    uniform unsigned int64 nc;
-    uniform unsigned int64 nb_panels;
+    uniform uint64 ns;
+    uniform uint64 nc;
+    uniform uint64 nb_panels;
     uniform SoA3D v; // vertices
     uniform SoA3D colloc; // collocation points
     uniform SoA3D normal; // normals
@@ -85,12 +85,12 @@ inline void kernel_symmetry(float3& inf, float3 colloc, const uniform float3& ve
 export void kernel_influence(
     uniform const MeshProxy& m,
     uniform float* uniform lhs,
-    uniform unsigned int64 ia, uniform unsigned int64 lidx, uniform float sigma
+    uniform uint64 ia, uniform uint64 lidx, uniform float sigma
     ) {
-    const uniform unsigned int64 v0 = lidx + lidx / m.ns;
-    const uniform unsigned int64 v1 = v0 + 1;
-    const uniform unsigned int64 v3 = v0 + m.ns+1;
-    const uniform unsigned int64 v2 = v3 + 1;
+    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]};
@@ -115,18 +115,18 @@ export uniform float kernel_trefftz_cd(
     uniform const MeshProxy& m,
     uniform float* uniform gamma,
     uniform float* uniform trefftz_buffer,
-    uniform unsigned int64 j, uniform unsigned int64 n, uniform float sigma
+    uniform uint64 j, uniform uint64 n, uniform float sigma
     ) {
-    uniform unsigned int64 begin = m.nb_panels + j;
-    uniform unsigned int64 end = begin + n;
+    uniform uint64 begin = m.nb_panels + j;
+    uniform uint64 end = begin + n;
     float cd = 0.0f;
 
     // Compute the induced velocity of the streamwise wake vortex segments
-    for (uniform unsigned int64 ia = m.nb_panels; ia < m.nb_panels + m.ns; ia++) {
-        const uniform unsigned int64 v0 = ia + ia / m.ns;
-        const uniform unsigned int64 v1 = v0 + 1;
-        const uniform unsigned int64 v3 = v0 + m.ns+1;
-        const uniform unsigned int64 v2 = v3 + 1;
+    for (uniform uint64 ia = m.nb_panels; ia < m.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 v2 = v3 + 1;
 
         // Broadcast vertices
         const uniform float3 vertex0 = {m.v.x[v0], m.v.y[v0], m.v.z[v0]};
@@ -147,7 +147,7 @@ export uniform float kernel_trefftz_cd(
     }
     // Perform the integration (Katz Plotkin, Low speed Aero | Eq 8.146)
     foreach(i = j ... j + n) {
-        unsigned int64 li = (m.nc-1) * m.ns + i;
+        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]};
         float dl = length(v1 - v0);
@@ -159,11 +159,11 @@ export uniform float kernel_trefftz_cd(
 export uniform float kernel_trefftz_cl(
     uniform const MeshProxy& m,
     uniform float* uniform gamma,
-    uniform unsigned int64 j, uniform unsigned int64 n
+    uniform uint64 j, uniform uint64 n
     ) {
     float cl = 0.0f;
     foreach(i = j ... j + n) {
-        unsigned int64 li = (m.nc-1) * m.ns + i;
+        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]};
         float dl = length(v1 - v0);

vlm/backends/cuda/src/vlm_backend_cuda.cu

@@ -134,6 +134,62 @@ void BackendCUDA::reset() {
     default_backend.reset();
 }
 
+struct SoA3D {
+    float* x;
+    float* y;
+    float* z;
+};
+
+struct MeshProxy {
+    uint64_t ns;
+    uint64_t nc;
+    uint64_t nb_panels;
+    SoA3D v; // vertices
+    SoA3D colloc; // collocation points
+    SoA3D normal; // normals
+};
+
+__global__ void kernel_influence_cuda(const MeshProxy& m, float* d_lhs, uint64_t lidx, float sigma) {
+
+}
+
+// void vlm::kernel_influence(
+//     u64 m, u64 n,
+//     f32 lhs[],
+//     f32 vx[], f32 vy[], f32 vz[],
+//     f32 collocx[], f32 collocy[], f32 collocz[],
+//     f32 normalx[], f32 normaly[], f32 normalz[],
+//     f32 sigma
+//     ) {
+
+//     // parallel for
+//     for (u64 lidx = 0; lidx < m*n; lidx++) {
+//         const u64 nb_panels = m * n;
+//         const u64 v0 = lidx + lidx / n;
+//         const u64 v1 = v0 + 1;
+//         const u64 v3 = v0 + n+1;
+//         const u64 v2 = v3 + 1;
+
+//         float3 vertex0{vx[v0], vy[v0], vz[v0]};
+//         float3 vertex1{vx[v1], vy[v1], vz[v1]};
+//         float3 vertex2{vx[v2], vy[v2], vz[v2]};
+//         float3 vertex3{vx[v3], vy[v3], vz[v3]};
+
+//         for (u64 ia2 = 0; ia2 < nb_panels; ia2++) {
+//             const float3 colloc(collocx[ia2], collocy[ia2], collocz[ia2]);
+//             float3 inf(0.0f, 0.0f, 0.0f);
+//             const float3 normal(normalx[ia2], normaly[ia2], normalz[ia2]);
+
+//             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);
+//             // store in col major order
+//             lhs[lidx * nb_panels + ia2] += linalg::dot(inf, normal);
+//         }
+//     }
+// }
+
 void BackendCUDA::compute_lhs(const FlowData& flow) {
     default_backend.compute_lhs(flow);
 }