From fe2e776076d88bfc53dcb6e525f0d48dd6fbf0fb Mon Sep 17 00:00:00 2001
From: Orestis <orestis.malaspinas@pm.me>
Date: Fri, 19 Feb 2021 21:23:00 +0100
Subject: [PATCH] forgot some commits
---
apps/kida/kida.c | 13 ++++++++++++-
meson.build | 3 ++-
src/c_interface/base64.c | 8 ++++----
src/c_interface/globals.h | 28 +++++++++++++++++++++++++++
src/c_interface/includes.h | 2 ++
src/c_interface/lbm_utils.c | 2 +-
src/c_interface/meson.build | 2 ++
src/c_interface/point.c | 4 ++--
src/c_interface/point.h | 4 ++--
src/c_interface/scalar_field.c | 35 ++++++++++++++++++++++++++++++++++
src/c_interface/scalar_field.h | 4 ++++
src/c_interface/triangle.c | 8 ++++++++
src/c_interface/triangle.h | 1 +
src/dynamics.fut | 6 ++++++
src/futhark.pkg | 1 +
src/liblbm_f32.fut | 3 +++
src/liblbm_f64.fut | 3 +++
src/meson.build | 31 +++++++++++++++++++++++++++---
18 files changed, 144 insertions(+), 14 deletions(-)
diff --git a/apps/kida/kida.c b/apps/kida/kida.c
index da31ead..d327e81 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 b5cafe1..c1ff24f 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 1c639d1..ef2fc15 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 5918918..0a93e90 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 48f471a..c22a711 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 a3194bc..6c6f703 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 8dbc4cb..266a697 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 b7ef7e5..8cbb916 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 c8e9946..bcd5b3a 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 656a135..2517ba2 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 29fe911..76375b2 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 8ce853f..5ef570c 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 063984c..29b1727 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 05b1fe3..c1c1854 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 c05fea5..4b152a1 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 0e7309b..09762b4 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 7f8f6c4..20d8830 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 ee65f6f..dbef9df 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'])
--
GitLab