diff --git a/apps/kida/kida.c b/apps/kida/kida.c
index da31ead55054e00bf4a6ba449908e87a5fce5844..d327e813e1c5141a8a1e78cbb1c58bb9e6ddda84 100644
--- a/apps/kida/kida.c
+++ b/apps/kida/kida.c
@@ -92,10 +92,20 @@ int main(int argc, char **argv) {
if (prefix != NULL) {
printf("output at iter %d\n", it);
tensor_field_from_futhark_f32_4d_inplace(f, ctx, lattice);
- tensor_field_from_futhark_f32_4d_inplace(d, ctx, distances);
+ // tensor_field_from_futhark_f32_4d_inplace(d, ctx, distances);
lbm_compute_density_inplace(lattice, density);
lbm_compute_velocity_inplace(lattice, velocity);
+ triangle_field tri = triangle_field_polygonise(density, ctx, 1.0);
+
+ char *stl_fname = stl_create_fname(prefix, "iso_surface", it);
+ stl_file stl = stl_file_create(stl_fname);
+ stl_file_write_triangle_field(&stl, tri);
+ stl_file_close(&stl);
+ free(stl_fname);
+ triangle_field_deallocate(tri);
+
+
char *fname = vtk_create_fname(prefix, "out", it);
box domain = box_create(0, lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1);
@@ -105,6 +115,7 @@ int main(int argc, char **argv) {
vtk_file_write_scalar_field(&vtk, density, domain, "density", base64);
vtk_file_write_tensor_field(&vtk, velocity, domain, "velocity", base64);
vtk_file_close(&vtk);
+
free(fname);
}
diff --git a/meson.build b/meson.build
index b5cafe1dbfb4ad5652110ebee4a6a613f9807c78..c1ff24f71ccbff799eb4dad3ba1c3cf0d992e05a 100644
--- a/meson.build
+++ b/meson.build
@@ -1,6 +1,6 @@
project('palathark', 'c',
version : '0.1',
- default_options : ['warning_level=3', 'buildtype=debugoptimized', 'b_ndebug=if-release'])
+ default_options : ['warning_level=3', 'buildtype=debug', 'b_sanitize=address,undefined', 'b_ndebug=if-release'])
lattice = get_option('lattice')
precision = get_option('precision')
@@ -18,6 +18,7 @@ else
liblbm = 'liblbm_f64.fut'
endif
+
math_lib = meson.get_compiler('c').find_library('m')
math_dep = declare_dependency(dependencies: [math_lib])
diff --git a/src/c_interface/base64.c b/src/c_interface/base64.c
index 1c639d16ad2be80956fbe47c70a3f94688445ab4..ef2fc157218cb638908e2a37449404ec83a68085 100644
--- a/src/c_interface/base64.c
+++ b/src/c_interface/base64.c
@@ -56,11 +56,11 @@ unsigned char * base64_encode(const unsigned char *src, size_t len,
/*
// Modification: These lines are commented so that paraview
// can read the output.
+ */
// if (line_len >= 72) {
// *pos++ = '\n';
// line_len = 0;
// }
- */
}
if (end - in) {
@@ -79,10 +79,10 @@ unsigned char * base64_encode(const unsigned char *src, size_t len,
/*
// Modification: These lines are commented so that paraview
- // can read the output.
+ // can read the output. */
// if (line_len)
- // *pos++ = '\n';
- */
+ // *pos++ = '\n';
+
*pos = '\0';
if (out_len)
diff --git a/src/c_interface/globals.h b/src/c_interface/globals.h
index 59189189b57d81ad3e786a5db3d9f890db1336f0..0a93e905c172f61f68627078d4005f02a4f52cc4 100644
--- a/src/c_interface/globals.h
+++ b/src/c_interface/globals.h
@@ -14,21 +14,49 @@ typedef double T;
static const T tol = DBL_EPSILON;
static const T min_val = DBL_MIN;
+
+typedef struct futhark_f64_2d futhark_2d;
+typedef struct futhark_f64_3d futhark_3d;
typedef struct futhark_f64_4d futhark_4d;
+
#define futhark_shape_4d(a, b) futhark_shape_f64_4d(a, b)
#define futhark_free_4d(a, b) futhark_free_f64_4d(a, b)
#define futhark_values_4d(a, b, c) futhark_values_f64_4d(a, b, c)
#define futhark_new_4d(a, b, c, d, e, f) futhark_new_f64_4d(a, b, c, d, e, f)
+
+#define futhark_shape_3d(a, b) futhark_shape_f64_3d(a, b)
+#define futhark_free_3d(a, b) futhark_free_f64_3d(a, b)
+#define futhark_values_3d(a, b, c) futhark_values_f64_3d(a, b, c)
+#define futhark_new_3d(a, b, c, d, e) futhark_new_f64_3d(a, b, c, d, e)
+
+#define futhark_shape_2d(a, b) futhark_shape_f64_2d(a, b)
+#define futhark_free_2d(a, b) futhark_free_f64_2d(a, b)
+#define futhark_values_2d(a, b, c) futhark_values_f64_2d(a, b, c)
+#define futhark_new_2d(a, b, c, d) futhark_new_f64_2d(a, b, c, d)
#else
typedef float T;
static const T tol = FLT_EPSILON;
static const T min_val = FLT_MIN;
+
+typedef struct futhark_f32_2d futhark_2d;
+typedef struct futhark_f32_3d futhark_3d;
typedef struct futhark_f32_4d futhark_4d;
+
#define futhark_shape_4d(a, b) futhark_shape_f32_4d(a, b)
#define futhark_free_4d(a, b) futhark_free_f32_4d(a, b)
#define futhark_values_4d(a, b, c) futhark_values_f32_4d(a, b, c)
#define futhark_new_4d(a, b, c, d, e, f) futhark_new_f32_4d(a, b, c, d, e, f)
+
+#define futhark_shape_3d(a, b) futhark_shape_f32_3d(a, b)
+#define futhark_free_3d(a, b) futhark_free_f32_3d(a, b)
+#define futhark_values_3d(a, b, c) futhark_values_f32_3d(a, b, c)
+#define futhark_new_3d(a, b, c, d, e) futhark_new_f32_3d(a, b, c, d, e)
+
+#define futhark_shape_2d(a, b) futhark_shape_f32_2d(a, b)
+#define futhark_free_2d(a, b) futhark_free_f32_2d(a, b)
+#define futhark_values_2d(a, b, c) futhark_values_f32_2d(a, b, c)
+#define futhark_new_2d(a, b, c, d) futhark_new_f32_2d(a, b, c, d)
#endif
static const T pi = M_PI;
diff --git a/src/c_interface/includes.h b/src/c_interface/includes.h
index 48f471aa26d3dd60b6883749c9cc29c3352f3b9d..c22a711024c1b0203b723ea3088722acbcf8f318 100644
--- a/src/c_interface/includes.h
+++ b/src/c_interface/includes.h
@@ -10,10 +10,12 @@
#include "point.h"
#include "box.h"
#include "vtk.h"
+#include "stl.h"
#include "units.h"
#include "utils.h"
#include "point.h"
#include "triangle.h"
+#include "triangle_field.h"
#include "base64.h"
#endif
\ No newline at end of file
diff --git a/src/c_interface/lbm_utils.c b/src/c_interface/lbm_utils.c
index a3194bc8de4f570486b185d39897a611b79fee30..6c6f7030d8d02fec1c685e1db03cb10ddbb5d397 100644
--- a/src/c_interface/lbm_utils.c
+++ b/src/c_interface/lbm_utils.c
@@ -17,7 +17,7 @@ tensor_field create_distances(fint nx, fint ny, fint nz, fint q, triangle *t, fi
for (fint iq = 0; iq < q; ++iq) {
T distance = (T)1;
point loc_point = point_create(ix, iy, iz);
- point next_vec = point_create_from_array_int(c_t[iq]);
+ point next_vec = point_from_array_int(c_t[iq]);
// printf("iq = %d, ",iq);
// point_print(next_vec);
// printf("\n");
diff --git a/src/c_interface/meson.build b/src/c_interface/meson.build
index 8dbc4cbfa8855cb36f080d63c5c6679156e2e1ab..266a697e7cd0dddf36a6c795a910820191ae401d 100644
--- a/src/c_interface/meson.build
+++ b/src/c_interface/meson.build
@@ -6,10 +6,12 @@ c_interface_sources = [
'units.c',
'tensor_field_2d.c',
'tensor_field.c',
+ 'triangle_field.c',
'box.c',
'triangle.c',
'point.c',
'vtk.c',
+ 'stl.c',
'scalar_field_2d.c']
c_interface = static_library('c_interface',
diff --git a/src/c_interface/point.c b/src/c_interface/point.c
index b7ef7e587a9c8b9633b459fb63d2b7a746aaece4..8cbb9161168ab3d6bd5f61b677a34305bba65b2f 100644
--- a/src/c_interface/point.c
+++ b/src/c_interface/point.c
@@ -16,7 +16,7 @@ point point_create_default() {
return point_create((T)0, (T)0, (T)0);
}
-point point_create_from_array_int(const fint array[3]) {
+point point_from_array_int(const fint array[3]) {
point p;
p.x = (T)array[0];
p.y = (T)array[1];
@@ -25,7 +25,7 @@ point point_create_from_array_int(const fint array[3]) {
return p;
}
-point point_create_from_array(const T array[3]) {
+point point_from_array(const T array[3]) {
point p;
p.x = array[0];
p.y = array[1];
diff --git a/src/c_interface/point.h b/src/c_interface/point.h
index c8e994640a7174d6b052abc187b56e41765d1d35..bcd5b3a590ca04f11a97ab4bf741337a7c381e7f 100644
--- a/src/c_interface/point.h
+++ b/src/c_interface/point.h
@@ -10,8 +10,8 @@ typedef struct {
point point_create(T x, T y, T z);
point point_create_default();
-point point_create_from_array_int(const fint array[3]);
-point point_create_from_array(const T array[3]);
+point point_from_array_int(const fint array[3]);
+point point_from_array(const T array[3]);
bool point_is_approx_equal(point lhs, point rhs, T tol);
point point_neg(point rhs);
diff --git a/src/c_interface/scalar_field.c b/src/c_interface/scalar_field.c
index 656a13537ee5a0b32ef2e1f8795c5471230b154c..2517ba2f6943a5ad673a1e03476368aa7ed2948e 100644
--- a/src/c_interface/scalar_field.c
+++ b/src/c_interface/scalar_field.c
@@ -146,6 +146,41 @@ void scalar_field_write_file(scalar_field field, char *dir, char *prefix, int i)
fclose(fp);
}
+void scalar_field_from_futhark_f32_3d_inplace(futhark_3d *f, struct futhark_context *ctx, scalar_field field) {
+#ifndef NDEBUG
+ const int64_t *shape = futhark_shape_3d(ctx, f);
+ fint nx = shape[0];
+ fint ny = shape[1];
+ fint nz = shape[2];
+#endif
+
+ assert(nx == field.nx);
+ assert(ny == field.ny);
+ assert(nz == field.nz);
+
+ futhark_values_3d(ctx, f, field.data);
+}
+
+scalar_field scalar_field_from_futhark_f32_3d(futhark_3d *f, struct futhark_context *ctx) {
+ const int64_t *shape = futhark_shape_3d(ctx, f);
+ fint nx = shape[0];
+ fint ny = shape[1];
+ fint nz = shape[2];
+
+ scalar_field field = scalar_field_allocate(nx, ny, nz);
+
+ scalar_field_from_futhark_f32_3d_inplace(f, ctx, field);
+
+ return field;
+}
+
+futhark_3d *scalar_field_to_futhark_f32_3d(scalar_field field, struct futhark_context *ctx) {
+
+ futhark_3d *f = futhark_new_3d(ctx, field.data, field.nx, field.ny, field.nz);
+ assert(f != NULL);
+ return f;
+}
+
// TODO: add possibility to have rho/u form futhark directly
// void scalar_field_from_futhark_f32_3d_inplace(struct futhark_f32_3d *f, struct futhark_context *ctx, scalar_field field) {
// assert(f != NULL);
diff --git a/src/c_interface/scalar_field.h b/src/c_interface/scalar_field.h
index 29fe91155644dec12b29234e0168ec971c69206a..76375b2064c093ea25664fd03bcd932359cfa55e 100644
--- a/src/c_interface/scalar_field.h
+++ b/src/c_interface/scalar_field.h
@@ -37,6 +37,10 @@ scalar_field_2d scalar_field_extract_2d_y(scalar_field field, int iy);
scalar_field_2d scalar_field_extract_2d_z(scalar_field field, int iz);
+void scalar_field_from_futhark_f32_3d_inplace(futhark_3d *rho, struct futhark_context *ctx, scalar_field field);
+scalar_field scalar_field_from_futhark_f32_3d(futhark_3d *rho, struct futhark_context *ctx);
+futhark_3d *scalar_field_to_futhark_f32_3d(scalar_field field, struct futhark_context *ctx);
+
// TODO: add possibility to have rho/u form futhark directly
// void scalar_field_from_futhark_f32_3d_inplace(struct futhark_f32_3d *f, struct futhark_context *ctx, scalar_field field);
diff --git a/src/c_interface/triangle.c b/src/c_interface/triangle.c
index 8ce853fbbcecc52c73f4e099b645e65d44999aa0..5ef570c1c2ddf3ac786f175560744801c7e1bf8c 100644
--- a/src/c_interface/triangle.c
+++ b/src/c_interface/triangle.c
@@ -133,6 +133,14 @@ void triangle_print(triangle t) {
point_print(t.p2);
}
+triangle triangle_from_array(const T data[9]) {
+ point p0 = point_from_array(data);
+ point p1 = point_from_array(&data[3]);
+ point p2 = point_from_array(&data[6]);
+
+ return triangle_create(p0, p1, p2);
+}
+
static char *stl_create_fname(char const * const prefix, char *name) {
char *buffer = malloc(80 * sizeof(char));
sprintf(buffer, "%s%s.stl", prefix, name);
diff --git a/src/c_interface/triangle.h b/src/c_interface/triangle.h
index 063984c4d2b6e7cbbc1950a4c9c5398081632228..29b1727164941a03e2c8264b77ded021fb06eec2 100644
--- a/src/c_interface/triangle.h
+++ b/src/c_interface/triangle.h
@@ -10,6 +10,7 @@ typedef struct {
} triangle;
triangle triangle_create(point p0, point p1, point p2);
+triangle triangle_from_array(const T data[9]);
bool triangle_is_approx_equal(triangle lhs, triangle rhs, T eps);
diff --git a/src/dynamics.fut b/src/dynamics.fut
index 05b1fe383e06a043ad287d88e35893f4363ea580..c1c1854098320ca8a0be52adcc9897c0f0d8a562 100644
--- a/src/dynamics.fut
+++ b/src/dynamics.fut
@@ -48,6 +48,9 @@ module mk_bgk (real: real) (d: descriptor with real = real.t) : (dynamics with r
let u = d.map_d( (real.f64 1) / rho * ) j
let feq = compute_feq rho u
in d.map2_q(\fi feqi -> fi * ( (real.f64 1) - omega) + feqi * omega) f feq
+ -- mad function used here. Not sure it adds any efficiency
+ -- in d.map2_q(\fi feqi -> real.mad (fi) ( (real.f64 1) - omega) (feqi * omega)) f feq
+
let compute_velocity (f: d.vec_q real) =
let j = compute_j f
@@ -164,6 +167,9 @@ module mk_reg (real: real)
let feq = compute_feq rho u
let fneq = compute_fneq f feq
let fneq = regularize fneq u
+
in d.map2_q(\feqi fneqi -> fneqi * ( (real.f64 1) - omega) + feqi) feq fneq
+ -- mad function used here. Not sure it adds any efficiency
+ -- in d.map2_q(\feqi fneqi -> real.mad (fneqi) ( (real.f64 1) - omega) (feqi)) feq fneq
}
diff --git a/src/futhark.pkg b/src/futhark.pkg
index c05fea53ef883b0070e53c8c6dd7e30dc0c9bf8c..4b152a13d60a62d51bbdca7f44ae638c8f7fa91b 100644
--- a/src/futhark.pkg
+++ b/src/futhark.pkg
@@ -1,3 +1,4 @@
require {
github.com/athas/vector 0.7.0 #f5aaa219a02949d1eb2978f81419cf5614abfa3d
+ github.com/omalaspinas/futcubes 0.0.0-20210219193616+d16367e9f96dbd7938d7916483386a8aa0516723 #d16367e9f96dbd7938d7916483386a8aa0516723
}
diff --git a/src/liblbm_f32.fut b/src/liblbm_f32.fut
index 0e7309bb533c72456a2677a2ccb15a9fe8785cce..09762b4df7dc4fc89ec87c9ccc4ea1699e7a0742 100644
--- a/src/liblbm_f32.fut
+++ b/src/liblbm_f32.fut
@@ -1,6 +1,7 @@
-- Program for Lattice-Boltzmann in Futhark
import "globals_f32"
import "sim"
+import "polygonise"
module sim_d3q27_reg_bfl_f32 = mk_sim d3q27_reg_bfl_f32
@@ -10,6 +11,8 @@ module sim_d3q19_reg_bfl_f32 = mk_sim d3q19_reg_bfl_f32
entry main_d3q19_reg = sim_d3q19_reg_bfl_f32.run
+entry main_polygonise = polygonise_f32
+
-- Test the sub1 function.
-- ==
-- input { 21i64 21i64 21i64 27i64 3i64 1.5f32 1000i64 }
diff --git a/src/liblbm_f64.fut b/src/liblbm_f64.fut
index 7f8f6c46ea49086e446f9463cf297633463277ee..20d8830fc1c2e13d48abce20233d7eb077889489 100644
--- a/src/liblbm_f64.fut
+++ b/src/liblbm_f64.fut
@@ -1,6 +1,7 @@
-- Program for Lattice-Boltzmann in Futhark
import "globals_f64"
import "sim"
+import "polygonise"
module sim_d3q27_reg_bfl_f64 = mk_sim d3q27_reg_bfl_f64
@@ -10,6 +11,8 @@ module sim_d3q19_reg_bfl_f64 = mk_sim d3q19_reg_bfl_f64
entry main_d3q19_reg_f64 = sim_d3q19_reg_bfl_f64.run
+entry main_polygonise = polygonise_f64
+
-- -- executes Lattice-Boltzmann it times and
-- let lattice_boltzmann [nx][ny][nz] (fin: [nx][ny][nz](lbm_model.vec_q real)) (vel: [nx][ny][nz](lbm_model.vec_d real)) (dists: [nx][ny][nz](lbm_model.vec_q real)) (omega: real) (max_t: i64): [nx][ny][nz](lbm_model.vec_q real) =
diff --git a/src/meson.build b/src/meson.build
index ee65f6f30f4273e82521640ef41a89e6f3b9f225..dbef9df9df115eb9ca7a036def38170247cfd8dc 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -1,15 +1,38 @@
+fut_sources = [
+ 'bc.fut',
+ 'vectors.fut',
+ 'utilities.fut',
+ 'sym_tens.fut',
+ 'sim.fut',
+ 'liblbm_f64.fut',
+ 'liblbm_f32.fut',
+ 'lbm.fut',
+ 'lattice.fut',
+ 'globals.fut',
+ 'globals_f64.fut',
+ 'globals_f32.fut',
+ 'descriptor.fut',
+ 'data_analysis.fut',
+ 'dynamics.fut',
+ 'polygonise.fut']
+
+fut_dep = declare_dependency(
+ sources : fut_sources
+)
+
futhark_generated = custom_target(
'futhark',
- input: join_paths(meson.source_root(),'src',liblbm),
+ input: [join_paths(meson.source_root(), 'src', liblbm), fut_sources],
output: ['liblbm.c', 'liblbm.h'],
- command: [futhark, futhark_backend, '@INPUT@', '--library', '-o', 'src/liblbm'],
+ command: [futhark, futhark_backend, '@INPUT0@', '--library', '-o', 'src/liblbm'],
install: true,
- install_dir: 'src'
+ install_dir: 'src',
)
liblbm_static = static_library('liblbm_static',
[futhark_generated],
c_args : ['-Wall', '-Wextra', '-pedantic', build_args],
+ include_directories: './',
)
liblbm_dep = declare_dependency(
@@ -17,6 +40,8 @@ liblbm_dep = declare_dependency(
link_with : [liblbm_static],
)
+
+
# gen = generator(futhark,
# output: ['@BASENAME@.c', '@BASENAME@.h'],
# arguments : [futhark_backend, '@INPUT@', '--library', '-o', 'src/liblbm'])