Added LBM from @orestis.malaspin

futhark_mpi/dispatch.c

@@ -1057,6 +1057,10 @@ static envelope_t *get_outer_envelope_2d(struct dispatch_context *dc, int thickn
     return outer_envelope;
 }
 
+typedef struct lbm_values {
+    float values[27];
+} lbm_values_t;
+
 static envelope_t *get_outer_envelope_3d(struct dispatch_context *dc, int thickness, envelope_t *inner_envelope) {
     int coordinate_y = dc->coordinates[0];
     int coordinate_x = dc->coordinates[1];
@@ -1990,17 +1994,17 @@ get_chunk_with_envelope(struct dispatch_context *dc, struct futhark_context *fc,
     switch (dc->n_dimensions) {
         case 1:
             result = get_chunk_with_envelope_1d(dc, fc, outer_envelope,
-                                              (void *(*)(struct futhark_context *, const void *, int64_t)) f);
+                                                (void *(*)(struct futhark_context *, const void *, int64_t)) f);
             break;
         case 2:
             result = get_chunk_with_envelope_2d(dc, fc, outer_envelope,
-                                              (void *(*)(struct futhark_context *, const void *, int64_t,
-                                                         int64_t)) f);
+                                                (void *(*)(struct futhark_context *, const void *, int64_t,
+                                                           int64_t)) f);
             break;
         case 3:
             result = get_chunk_with_envelope_3d(dc, fc, outer_envelope,
-                                              (void *(*)(struct futhark_context *, const void *, int64_t, int64_t,
-                                                         int64_t)) f);
+                                                (void *(*)(struct futhark_context *, const void *, int64_t, int64_t,
+                                                           int64_t)) f);
             break;
     }
     envelope_free(outer_envelope);

futhark_mpi/dispatch.h

@@ -60,10 +60,6 @@ extern envelope_t *get_inner_envelope(struct dispatch_context *dc, struct futhar
 
 extern envelope_t *get_outer_envelope(struct dispatch_context *dc, struct futhark_context *fc, int thickness);
 
-//extern void *
-//futhark_outer_envelope_new(struct dispatch_context *dc, struct futhark_context *fc, envelope_t *outer_envelope,
-//                           void *f(struct futhark_context *, const void *), char *futhark_type);
-
 extern chunk_info_t get_chunk_info(struct dispatch_context *dc);
 
 extern void *get_chunk_with_envelope(struct dispatch_context *dc, struct futhark_context *fc, int thickness,

futhark_mpi/envelope.fut

@@ -163,8 +163,8 @@ let augment_2d [n][m][ty][tx] (matrix:[n][m]u8) (envelope: envelope_2d_t [n][m][
     ) :> [new_n][new_m]u8
 
 let augment_3d [n][m][l] (cube :[n][m][l]u8) (envelope: envelope_3d_t) :[][][]u8 =
-    let ty = length envelope.top.south_west
-    let tx = length envelope.top.south_west[0]
+    let ty = length envelope.top.surface
+    let tx = length envelope.left.surface[0]
     let tz = length envelope.front.surface[0,0]
     let new_n = n + ty * 2
     let new_m = m + tx * 2
@@ -179,7 +179,7 @@ let augment_3d [n][m][l] (cube :[n][m][l]u8) (envelope: envelope_3d_t) :[][][]u8
         -- Top-South
         else if (i >= 0 && i < ty && j >= tx && j < (m+tx) && k >= 0 && k < tz) then envelope.top.south[i,j-tx,k]
         -- Top-Surface
-        else if (i >= 0 && i < ty && j >= tx && j < (m+tx) && k >= tz && k < (l+tz)) then envelope.top.surface[i,j-tx, k-tz]
+        else if (i >= 0 && i < ty && j >= tx && j < (m+tx) && k >= tz && k < (l+tz)) then envelope.top.surface[i,j-tx,k-tz]
         -- Top North
         else if (i >= 0 && i < ty && j >= tx && j < (m+tx) && k >= (l+tz) && k < new_l) then envelope.top.north[i,j-tx, k-l-tz]
         -- Top South East
@@ -187,7 +187,7 @@ let augment_3d [n][m][l] (cube :[n][m][l]u8) (envelope: envelope_3d_t) :[][][]u8
         -- Top East
         else if (i >= 0 && i < ty && j >= m+tx && j < new_m && k >= tz && k < l + tz) then envelope.top.east[i,j-m-tx,k-tz]
         -- Top North East
-        else if (i >= 0 && i < ty && j >= m+tx && j < new_m && k >= k+tz && k < new_l) then envelope.top.north_east[i,j-m-tx,k-l-tz]
+        else if (i >= 0 && i < ty && j >= m+tx && j < new_m && k >= l+tz && k < new_l) then envelope.top.north_east[i,j-m-tx,k-l-tz]
 
         -- Front-West
         else if (i >= ty && i < n+ty && j >= 0 && j < tx && k >= 0 && k < tz) then envelope.front.west[i-ty, j, k]
@@ -195,12 +195,16 @@ let augment_3d [n][m][l] (cube :[n][m][l]u8) (envelope: envelope_3d_t) :[][][]u8
         else if (i >= ty && i < n+ty && j >= 0 && j < tx && k >= tz && k < l+tz) then envelope.left.surface[i-ty, j, k-tz]
         -- Back West
         else if (i >= ty && i < n+ty && j >= 0 && j < tx && k >= l+tz && k < new_l) then envelope.back.west[i-ty, j, k-l-tz]
+        -- Front Surface
+        else if (i >= ty && i < n+ty && j >= tx && j < m+tx && k >= 0 && k < tz) then envelope.front.surface[i-ty, j-tx, k]
+        -- Back Surface
+        else if (i >= ty && i < n+ty && j >= tx && j < m+tx && k >= l+tz && k < new_l) then envelope.back.surface[i-ty, j-tx, k-l-tz]
         -- Front East
         else if (i >= ty && i < n+ty && j >= m+tx && j < new_m && k >= 0 && k < tz) then envelope.front.east[i-ty, j-m-tx, k]
         -- Right
-        else if (i >= ty && i < n+ty && j >= m+tx && j < new_m && k >= tz && k < l+tz) then envelope.left.surface[i-ty, j-m-tx, k-tz]
+        else if (i >= ty && i < n+ty && j >= m+tx && j < new_m && k >= tz && k < l+tz) then envelope.right.surface[i-ty, j-m-tx, k-tz]
         -- Back East
-        else if (i >= ty && i < n+ty && j >= m+tx && j < new_m && k >= l+tz && k < new_l) then envelope.back.west[i-ty, j-m-tx, k-l-tz]
+        else if (i >= ty && i < n+ty && j >= m+tx && j < new_m && k >= l+tz && k < new_l) then envelope.back.east[i-ty, j-m-tx, k-l-tz]
 
         -- Bottom-South-West
         else if (i >= n+ty && i < new_n && j >= 0 && j < tx && k >= 0 && k < tz) then envelope.bottom.south_west[i-n-ty,j,k]
@@ -219,7 +223,7 @@ let augment_3d [n][m][l] (cube :[n][m][l]u8) (envelope: envelope_3d_t) :[][][]u8
         -- Bottom East
         else if (i >= n+ty && i < new_n && j >= m+tx && j < new_m && k >= tz && k < l + tz) then envelope.bottom.east[i-n-ty,j-m-tx,k-tz]
         -- Bottom North East
-        else if (i >= n+ty && i < new_n && j >= m+tx && j < new_m && k >= k+tz && k < new_l) then envelope.bottom.north_east[i-n-ty,j-m-tx,k-l-tz]
+        else if (i >= n+ty && i < new_n && j >= m+tx && j < new_m && k >= l+tz && k < new_l) then envelope.bottom.north_east[i-n-ty,j-m-tx,k-l-tz]
 
-        else cube[i-ty, j-tx, l-tz]
+        else cube[i-ty, j-tx, k-tz]
     ) :> [new_n][new_m][new_l]u8

lattice_boltzmann/CMakeLists.txt

@@ -4,7 +4,7 @@ project(lattice_boltzmann C)
 set(CMAKE_C_STANDARD 11)
 
 if (CMAKE_BUILD_TYPE MATCHES Debug)
-    set(GCC_COMPILE_FLAGS "-Wall -Wextra -Wconversion -pedantic")
+    set(GCC_COMPILE_FLAGS "-Wall -Wextra -Wconversion -pedantic -fsanitize=address,undefined")
 elseif (CMAKE_BUILD_TYPE MATCHES Release)
     set(GCC_COMPILE_FLAGS "-O3")
 endif ()

lattice_boltzmann/lbm.fut

@@ -18,7 +18,117 @@ entry augment_1d [n][tx] (xs: [n]u8) (envelope: envelope_1d_t [tx]) :[]u8 = env.
 entry augment_2d [n][m][ty][tx] (matrix: [n][m]u8) (envelope: envelope_2d_t [n][m][tx][ty]) :[][]u8 = env.augment_2d matrix envelope
 entry augment_3d [n][m][l] (cube: [n][m][l]u8) (envelope: envelope_3d_t) :[][][]u8 = env.augment_3d cube envelope
 
+-- Program for Lattice-Boltzmann in Futhark
+-- Arrays everywhere
+
+-- Two useful functions for later
+let tabulate_4d 'a  (p: i64) (n: i64) (m: i64) (o: i64) (f: i64 -> i64 -> i64 -> i64 -> a): *[p][n][m][o]a =
+  map1 (f >-> tabulate_3d n m o) (iota p)
+
+let dotprod (xs: []f32) (ys: []f32): f32 =
+  reduce (+) 0 (map (\(x, y) -> x * y) (zip xs ys))
+
+let map3d 'a [nx][ny][nz]'b (f: a -> b) (xs: [nx][ny][nz]a): [nx][ny][nz]b =
+  map(\xs1 ->
+    map(\xs2 ->
+      map(\xs3 ->
+        f xs3
+      ) xs2
+    ) xs1
+  ) xs
+
+let map2_3d 'a 'c [nx][ny][nz]'b (f: a -> c -> b) (xs: [nx][ny][nz]a) (ys: [nx][ny][nz]c): [nx][ny][nz]b =
+  map2(\xs1 ys1 ->
+    map2(\xs2 ys2 ->
+      map2(\xs3 ys3 ->
+        f xs3 ys3
+      ) xs2 ys2
+    ) xs1 ys1
+  ) xs ys
+
+------------------------------------
+------------ CONSTANTS -------------
+------------------------------------
+
+module d3q27 = {
+
+    let c: [][]f32 = [
+                [ 0, 0, 0],
+
+                [-1, 0, 0], [ 0,-1, 0], [ 0, 0,-1],
+                [-1,-1, 0], [-1, 1, 0], [-1, 0,-1],
+                [-1, 0, 1], [ 0,-1,-1], [ 0,-1, 1],
+                [-1,-1,-1], [-1,-1, 1], [-1, 1,-1], [-1, 1, 1],
+
+                [ 1, 0, 0], [ 0, 1, 0], [ 0, 0, 1],
+                [ 1, 1, 0], [ 1,-1, 0], [ 1, 0, 1],
+                [ 1, 0,-1], [ 0, 1, 1], [ 0, 1,-1],
+                [ 1, 1, 1], [ 1, 1,-1], [ 1,-1, 1], [ 1,-1,-1] ]
+
+    let w: []f32 = [
+                8.0/27,
+                2.0/27,  2.0/27,  2.0/27,
+                1.0/54,  1.0/54,  1.0/54,
+                1.0/54,  1.0/54,  1.0/54,
+                1.0/216, 1.0/216, 1.0/216, 1.0/216,
+                2.0/27,  2.0/27,  2.0/27,
+                1.0/54,  1.0/54,  1.0/54,
+                1.0/54,  1.0/54,  1.0/54,
+                1.0/216, 1.0/216, 1.0/216, 1.0/216
+            ]
+}
+
+let compute_rho [nx][ny][nz][q] (f: [nx][ny][nz][q]f32): [nx][ny][nz]f32 =
+  map3d(\fxyz ->
+    reduce (+) 0 fxyz
+  ) f
+
+let compute_j [nx][ny][nz][q][d] (f: [nx][ny][nz][q]f32)(c: [q][d]f32): [nx][ny][nz][d]f32 =
+  map3d(\fxyz ->
+    map(\ci ->
+      dotprod ci fxyz
+    ) (transpose c) -- intrinsic
+  ) f
+
+let compute_feq [nx][ny][nz][q][d] (w: [q]f32) (c: [q][d]f32)(rho: [nx][ny][nz]f32)(j: [nx][ny][nz][d]f32): [nx][ny][nz][q]f32 =
+  map2_3d(\rho_xyz j_xyz ->
+    let u = map(\ji -> ji / rho_xyz) j_xyz
+    let u_sqr = dotprod u u
+
+    in map2(\wi ci ->
+      let c_u = dotprod ci u
+      in wi * rho_xyz* (1 + 3*c_u + 4.5 * c_u * c_u - 1.5 * u_sqr)
+    ) w c
+  ) rho j
+
+let collision [nx][ny][nz][q] (f: [nx][ny][nz][q]f32)(feq: [nx][ny][nz][q]f32) (omega: f32): [nx][ny][nz][q]f32 =
+  map2_3d (\f_xyz feq_xyz ->
+    map2(\f_i feq_i ->
+      (1-omega) * f_i + omega * feq_i
+    ) f_xyz feq_xyz
+  ) f feq
+
+let collide [nx][ny][nz][q][d] (f: [nx][ny][nz][q]f32) (omega: f32) (w: [q]f32) (c: [q][d]f32): [nx][ny][nz][q]f32 =
+  let rho = compute_rho f
+  let j = compute_j f c
+  let feq = compute_feq w c rho j
+
+  in collision f feq omega
+
+let streaming [nx][ny][nz][q][d] (f: [nx][ny][nz][q]f32) (c: [q][d]f32): [nx][ny][nz][q]f32 =
+    tabulate_4d nx ny nz q (\x y z i ->
+      let next_x = (x - (i64.f32 c[i,0]) + nx ) % nx
+      let next_y = (y - (i64.f32 c[i,1]) + ny ) % ny
+      let next_z = (z - (i64.f32 c[i,2]) + nz ) % nz
+
+      in f[next_x, next_y, next_z, i]
+    )
+
+let main_loop [nx][ny][nz] (fin: [nx][ny][nz][27]f32)(omega: f32) (max_t: i64): [nx][ny][nz][27]f32 =
+  let f = loop f = fin for _i < max_t do
+    streaming (collide f omega (copy d3q27.w) (copy d3q27.c)) (copy d3q27.c)
+  in f
+
 entry next_chunk_lbm [n][m][l][b] (chunk_lbm :[n][m][l][b]f32) :[][][][b]f32 =
-    -- TODO: Here, compute next lbm
-    let next_lbm = chunk_lbm
+    let next_lbm = (main_loop (chunk_lbm :> [n][m][l][27]f32) 1.0 1) :> [][][][b]f32
     in next_lbm[1:n-1, 1:m-1, 1:l-1, :]

lattice_boltzmann/main.c

@@ -6,6 +6,7 @@
  */
 
 #include <stdio.h>
+#include <stddef.h>
 #include "stdlib.h"
 #include <mpi.h>
 #include "lbm.h"
@@ -20,20 +21,20 @@ typedef struct lbm_values {
     float values[NB_VALUES];
 } lbm_values_t;
 
-
 void init_chunk_lbm(chunk_info_t *ci) {
-    float *data32 = ci->data;
+    lbm_values_t *data32 = ci->data;
     for (int i = 0; i < ci->dimensions[0]; ++i) {
         for (int j = 0; j < ci->dimensions[1]; ++j) {
             for (int k = 0; k < ci->dimensions[2]; ++k) {
                 int idx = INDEX_3D_TO_1D(i, j, k, ci->dimensions[1], ci->dimensions[2]);
-                data32[idx] = rand() % 2;
+                for (int l = 0; l < NB_VALUES; ++l) {
+                    data32[idx].values[l] = ((float) rand() / (float) rand());
+                }
             }
         }
     }
 }
 
-
 struct futhark_f32_4d *
 convert_chunk_with_envelope(struct futhark_context *fc, const void *data, int64_t dim0, int64_t dim1, int64_t dim2) {
     return futhark_new_f32_4d(fc, data, dim0, dim1, dim2, NB_VALUES);
@@ -70,21 +71,21 @@ int main(int argc, char *argv[]) {
 
     int lbm_dimensions[3] = {atoi(argv[1]), atoi(argv[2]), atoi(argv[3])};
 
-    MPI_Datatype cell;
-    int lengths[1] = {NB_VALUES};
-    MPI_Aint displacements[3];
-    struct lbm_values dummy_values;
-    MPI_Aint base_address;
-    MPI_Get_address(&dummy_values, &base_address);
-    displacements[0] = MPI_Aint_diff(displacements[0], base_address);
-
-    MPI_Datatype types[1] = {MPI_FLOAT};
-    MPI_Type_create_struct(1, lengths, displacements, types, &cell);
-    MPI_Type_commit(&cell);
-
-
-    struct dispatch_context *disp_context = dispatch_context_new(lbm_dimensions, cell, 3);
-    dispatch_context_print(disp_context);
+    // https://stackoverflow.com/a/33624425
+    int count = 1;
+    int array_of_blocklengths[] = {NB_VALUES};
+    MPI_Aint array_of_displacements[] = {offsetof(struct lbm_values, values)};
+    MPI_Datatype array_of_types[] = {MPI_FLOAT};
+    MPI_Datatype tmp_type, my_mpi_type;
+    MPI_Aint lb, extent;
+
+    MPI_Type_create_struct(count, array_of_blocklengths, array_of_displacements,
+                           array_of_types, &tmp_type);
+    MPI_Type_get_extent(tmp_type, &lb, &extent);
+    MPI_Type_create_resized(tmp_type, lb, extent, &my_mpi_type);
+    MPI_Type_commit(&my_mpi_type);
+
+    struct dispatch_context *disp_context = dispatch_context_new(lbm_dimensions, my_mpi_type, 3);
     chunk_info_t ci = get_chunk_info(disp_context);
     init_chunk_lbm(&ci);
 
@@ -95,5 +96,7 @@ int main(int argc, char *argv[]) {
     dispatch_context_free(disp_context);
     futhark_context_config_free(fut_config);
     futhark_context_free(fut_context);
+    MPI_Type_free(&tmp_type);
+    MPI_Type_free(&my_mpi_type);
     return MPI_Finalize();
 }