diff --git a/apps/cavity3d/cavity3d.c b/apps/cavity3d/cavity3d.c index adda2853fdc0ee487e1c95a1c342245e13f71f8a..a201e0fedb46ed9d74a023b8607622b6305480d4 100644 --- a/apps/cavity3d/cavity3d.c +++ b/apps/cavity3d/cavity3d.c @@ -1,12 +1,12 @@ -#include <liblbm.h> #include "c_interface/includes.h" +#include <assert.h> +#include <float.h> +#include <liblbm.h> +#include <math.h> #include <stdio.h> #include <stdlib.h> #include <string.h> -#include <assert.h> #include <time.h> -#include <float.h> -#include <math.h> #ifdef D3Q19 const fint l_q = 19; @@ -19,11 +19,13 @@ const fint l_q = 27; * Print how the program is supposed to be called */ void usage() { - printf("%s\n%s\n%s\n%s\n%s\n%s\n", "\nusage:\n ./lbm <n> <Re> <max_t> [<prefix>]", - "where:", " <n> is the mesh dimension in each spatial direction ( >= 10),", - "\t<Re> is the Reynolds number", - "\t<max_t> is maximum simulation time in physical units.", - "\t[<prefix>]: (optional argument) is the directory name where the output will be written."); + printf("%s\n%s\n%s\n%s\n%s\n%s\n", + "\nusage:\n ./lbm <n> <Re> <max_t> [<prefix>]", "where:", + " <n> is the mesh dimension in each spatial direction ( >= 10),", + "\t<Re> is the Reynolds number", + "\t<max_t> is maximum simulation time in physical units.", + "\t[<prefix>]: (optional argument) is the directory name where the " + "output will be written."); } int main(int argc, char **argv) { @@ -33,56 +35,66 @@ int main(int argc, char **argv) { return EXIT_FAILURE; } struct futhark_context_config *cfg = futhark_context_config_new(); - struct futhark_context *ctx = futhark_context_new(cfg); + struct futhark_context *ctx = futhark_context_new(cfg); printf("We are in the d = 3, q = %d\n", l_q); - const fint n = atoi(argv[1]); - const T Re = atof(argv[2]); - const T end_time = atof(argv[3]); + const fint n = atoi(argv[1]); + const T Re = atof(argv[2]); + const T end_time = atof(argv[3]); const char *prefix = NULL; if (5 == argc) { prefix = argv[4]; check_dir(prefix); } - const T log_t = 1.0/5.0; + const T log_t = 1.0 / 5.0; // const T log_t = 2.0; const T u_max = 0.1; - units_converter units = units_converter_create_inertial(Re, 1.0, n, 1.0, u_max, 1.0, 1.0, 1.0); + units_converter units = + units_converter_create_inertial(Re, 1.0, n, 1.0, u_max, 1.0, 1.0, 1.0); units_converter_print(units); const fint iter = (fint)(end_time / units.delta_t + 0.5); // Norme des vitesses en chaque point du maillage (struct futhark et c) T rho = 1.0; - // tensor_field lattice = lbm_allocate_and_ini_at_equilibrium_kida(rho, u_max, n+1, n+1, n+1, l_q); - tensor_field lattice = lbm_allocate_and_ini_at_equilibrium_cavity(rho, u_max, n+1, n+1, n+1, l_q); - - triangle *cavity = triangle_create_cube(point_create(0.5, 0.5, 0.5), point_create(n-0.5, n-0.5, n-0.5)); - tensor_field distances = create_distances(lattice.nx, lattice.ny, lattice.nz, lattice.q, cavity, 12); + // tensor_field lattice = lbm_allocate_and_ini_at_equilibrium_kida(rho, + // u_max, n+1, n+1, n+1, l_q); + tensor_field lattice = lbm_allocate_and_ini_at_equilibrium_cavity( + rho, u_max, n + 1, n + 1, n + 1, l_q); + + triangle *cavity = triangle_create_cube( + point_create(0.5, 0.5, 0.5), point_create(n - 0.5, n - 0.5, n - 0.5)); + tensor_field distances = create_distances( + lattice.nx, lattice.ny, lattice.nz, lattice.q, cavity, 12); if (prefix != NULL) { triangle_write_stl(cavity, 12, prefix, "cavity"); printf("output distances %d\n", 0); - char *fname = vtk_create_fname(prefix, "dist", 0); - vtk_file vtk = vtk_file_create_default(fname, box_create(0, lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1), 1.0); - vtk_file_write_tensor_field(&vtk, &distances, box_create(0, lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1), "distances", base64); + char *fname = vtk_create_fname(prefix, "dist", 0); + vtk_file vtk = vtk_file_create_default(fname, + box_create(0, lattice.nx - 1, 0, lattice.ny - 1, 0, lattice.nz - 1), + 1.0); + vtk_file_write_tensor_field(&vtk, &distances, + box_create(0, lattice.nx - 1, 0, lattice.ny - 1, 0, lattice.nz - 1), + "distances", base64); vtk_file_close(&vtk); - + free(fname); } - scalar_field density = lbm_compute_density(lattice); + scalar_field density = lbm_compute_density(lattice); tensor_field velocity = lbm_compute_velocity(lattice); - futhark_4d *f = tensor_field_to_futhark_f32_4d(lattice, ctx); - futhark_4d *d = tensor_field_to_futhark_f32_4d(distances, ctx); + futhark_4d *f = tensor_field_to_futhark_f32_4d(lattice, ctx); + futhark_4d *d = tensor_field_to_futhark_f32_4d(distances, ctx); futhark_4d *vel_ini = tensor_field_to_futhark_f32_4d(velocity, ctx); // TODO add version with rho/u as outputu from futhark // futhark_4d *f_vel = tensor_field_to_futhark_f32_4d(velocity, ctx); - // struct futhark_f32_3d *f_rho = scalar_field_to_futhark_f32_3d(density, ctx); + // struct futhark_f32_3d *f_rho = scalar_field_to_futhark_f32_3d(density, + // ctx); struct timespec start, finish; fint delta_iter = max((fint)(log_t / units.delta_t + 0.5), 1); @@ -101,32 +113,39 @@ int main(int argc, char **argv) { char *fname = vtk_create_fname(prefix, "out", it); - box domain = box_create(1, lattice.nx-2, 1, lattice.ny-2, 1, lattice.nz-2); + box domain = box_create( + 1, lattice.nx - 2, 1, lattice.ny - 2, 1, lattice.nz - 2); - // vtk_file vtk = vtk_file_create_default(fname, box_create(0, lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1), 1.0); + // vtk_file vtk = vtk_file_create_default(fname, box_create(0, + // lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1), 1.0); vtk_file vtk = vtk_file_create_default(fname, domain, 1.0); - vtk_file_write_scalar_field(&vtk, &density, domain, "density", base64); - vtk_file_write_tensor_field(&vtk, &velocity, domain, "velocity", base64); + 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); } - f = lbm_collide_and_stream(ctx, &f, vel_ini, d, units.omega, delta_iter); + f = lbm_collide_and_stream( + ctx, &f, vel_ini, d, units.omega, delta_iter); // TODO add version with rho/u as output - // f = lbm_collide_and_stream_rho_u(ctx, &f, &f_rho, &f_vel, units.omega, delta_iter); + // f = lbm_collide_and_stream_rho_u(ctx, &f, &f_rho, &f_vel, + // units.omega, delta_iter); i += 1; } - + futhark_context_sync(ctx); clock_gettime(CLOCK_MONOTONIC, &finish); T elapsed = (finish.tv_sec - start.tv_sec) * 1e6; elapsed += (finish.tv_nsec - start.tv_nsec) / 1e3; printf("elapsed microseconds = %f\n", elapsed); - printf("approximative MLUSP = %f\n", (T)units.nx * (T)units.ny * (T)units.nz * (T)i * (T)delta_iter / elapsed); + printf("approximative MLUSP = %f\n", (T)units.nx * (T)units.ny * + (T)units.nz * (T)i * + (T)delta_iter / elapsed); - if (prefix != NULL) { printf("output at iter %d\n", iter); tensor_field_from_futhark_f32_4d_inplace(f, ctx, lattice); @@ -135,18 +154,19 @@ int main(int argc, char **argv) { char *fname = vtk_create_fname(prefix, "out", iter); - box domain = box_create(1, lattice.nx-2, 1, lattice.ny-2, 1, lattice.nz-2); + box domain = + box_create(1, lattice.nx - 2, 1, lattice.ny - 2, 1, lattice.nz - 2); vtk_file vtk = vtk_file_create_default(fname, domain, 1.0); vtk_file_write_scalar_field(&vtk, &density, domain, "density", base64); - vtk_file_write_tensor_field(&vtk, &velocity, domain, "velocity", base64); + vtk_file_write_tensor_field( + &vtk, &velocity, domain, "velocity", base64); vtk_file_close(&vtk); - + free(fname); } free(cavity); - tensor_field_deallocate(lattice); tensor_field_deallocate(velocity); scalar_field_deallocate(density); @@ -155,12 +175,12 @@ int main(int argc, char **argv) { futhark_free_4d(ctx, f); futhark_free_4d(ctx, d); futhark_free_4d(ctx, vel_ini); -// TODO add version with rho/u output of futhark + // TODO add version with rho/u output of futhark // futhark_free_4d(ctx, f_vel); // futhark_free_f32_3d(ctx, f_rho); - futhark_context_config_free(cfg); futhark_context_free(ctx); + futhark_context_config_free(cfg); return EXIT_SUCCESS; } diff --git a/apps/kida/kida.c b/apps/kida/kida.c index 3dd736a4cd07593dbde3e464fcde9d0bfc4cf3bb..3dcaa28ea3e9c1ae613c080eeca370fbdcf8b6b0 100644 --- a/apps/kida/kida.c +++ b/apps/kida/kida.c @@ -1,12 +1,12 @@ -#include <liblbm.h> #include "c_interface/includes.h" +#include <assert.h> +#include <float.h> +#include <liblbm.h> +#include <math.h> #include <stdio.h> #include <stdlib.h> #include <string.h> -#include <assert.h> #include <time.h> -#include <float.h> -#include <math.h> #ifdef D3Q19 const fint l_q = 19; @@ -19,11 +19,13 @@ const fint l_q = 27; * Print how the program is supposed to be called */ void usage() { - printf("%s\n%s\n%s\n%s\n%s\n%s\n", "\nusage:\n ./lbm <n> <Re> <max_t> [<prefix>]", - "where:", " <n> is the mesh dimension in each spatial direction ( >= 10),", - "\t<Re> is the Reynolds number", - "\t<max_t> is maximum simulation time in physical units.", - "\t[<prefix>]: (optional argument) is the directory name where the output will be written."); + printf("%s\n%s\n%s\n%s\n%s\n%s\n", + "\nusage:\n ./lbm <n> <Re> <max_t> [<prefix>]", "where:", + " <n> is the mesh dimension in each spatial direction ( >= 10),", + "\t<Re> is the Reynolds number", + "\t<max_t> is maximum simulation time in physical units.", + "\t[<prefix>]: (optional argument) is the directory name where the " + "output will be written."); } int main(int argc, char **argv) { @@ -33,55 +35,63 @@ int main(int argc, char **argv) { return EXIT_FAILURE; } struct futhark_context_config *cfg = futhark_context_config_new(); - struct futhark_context *ctx = futhark_context_new(cfg); + struct futhark_context *ctx = futhark_context_new(cfg); futhark_context_config_set_debugging(cfg, 1); futhark_context_config_set_logging(cfg, 1); printf("We are in the d = 3, q = %d\n", l_q); - const fint n = atoi(argv[1]); - const T Re = atof(argv[2]); - const T end_time = atof(argv[3]); + const fint n = atoi(argv[1]); + const T Re = atof(argv[2]); + const T end_time = atof(argv[3]); const char *prefix = NULL; if (5 == argc) { prefix = argv[4]; check_dir(prefix); } - const T log_t = 1.0/5.0; + const T log_t = 1.0 / 5.0; // const T log_t = 2.0; const T u_max = 0.1; - units_converter units = units_converter_create_inertial(Re, 1.0, n, 1.0, u_max, 1.0, 1.0, 1.0); + units_converter units = + units_converter_create_inertial(Re, 1.0, n, 1.0, u_max, 1.0, 1.0, 1.0); units_converter_print(units); const fint iter = (fint)(end_time / units.delta_t + 0.5); // Norme des vitesses en chaque point du maillage (struct futhark et c) - T rho = 1.0; - tensor_field lattice = lbm_allocate_and_ini_at_equilibrium_kida(rho, u_max, n+1, n+1, n+1, l_q); + T rho = 1.0; + tensor_field lattice = lbm_allocate_and_ini_at_equilibrium_kida( + rho, u_max, n + 1, n + 1, n + 1, l_q); - tensor_field distances = create_distances_inactive(lattice.nx, lattice.ny, lattice.nz, lattice.q); + tensor_field distances = create_distances_inactive( + lattice.nx, lattice.ny, lattice.nz, lattice.q); if (prefix != NULL) { printf("output distances %d\n", 0); - char *fname = vtk_create_fname(prefix, "dist", 0); - vtk_file vtk = vtk_file_create_default(fname, box_create(0, lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1), 1.0); - vtk_file_write_tensor_field(&vtk, &distances, box_create(0, lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1), "distances", base64); + char *fname = vtk_create_fname(prefix, "dist", 0); + vtk_file vtk = vtk_file_create_default(fname, + box_create(0, lattice.nx - 1, 0, lattice.ny - 1, 0, lattice.nz - 1), + 1.0); + vtk_file_write_tensor_field(&vtk, &distances, + box_create(0, lattice.nx - 1, 0, lattice.ny - 1, 0, lattice.nz - 1), + "distances", base64); vtk_file_close(&vtk); - + free(fname); } - scalar_field density = lbm_compute_density(lattice); + scalar_field density = lbm_compute_density(lattice); tensor_field velocity = lbm_compute_velocity(lattice); - futhark_4d *f = tensor_field_to_futhark_f32_4d(lattice, ctx); - futhark_4d *d = tensor_field_to_futhark_f32_4d(distances, ctx); + futhark_4d *f = tensor_field_to_futhark_f32_4d(lattice, ctx); + futhark_4d *d = tensor_field_to_futhark_f32_4d(distances, ctx); futhark_4d *vel_ini = tensor_field_to_futhark_f32_4d(velocity, ctx); // TODO add version with rho/u as outputu from futhark // futhark_4d *f_vel = tensor_field_to_futhark_f32_4d(velocity, ctx); - // struct futhark_f32_3d *f_rho = scalar_field_to_futhark_f32_3d(density, ctx); + // struct futhark_f32_3d *f_rho = scalar_field_to_futhark_f32_3d(density, + // ctx); struct timespec start, finish; fint delta_iter = max((fint)(log_t / units.delta_t + 0.5), 1); @@ -100,50 +110,59 @@ int main(int argc, char **argv) { char *fname = vtk_create_fname(prefix, "out", it); - box domain = box_create(0, lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1); + box domain = box_create( + 0, lattice.nx - 1, 0, lattice.ny - 1, 0, lattice.nz - 1); - // vtk_file vtk = vtk_file_create_default(fname, box_create(0, lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1), 1.0); + // vtk_file vtk = vtk_file_create_default(fname, box_create(0, + // lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1), 1.0); vtk_file vtk = vtk_file_create_default(fname, domain, 1.0); - vtk_file_write_scalar_field(&vtk, &density, domain, "density", base64); - vtk_file_write_tensor_field(&vtk, &velocity, domain, "velocity", base64); + 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); - triangle_field tri = triangle_field_polygonise(density, ctx, 1.0, mtetras); // mtetras or mcubes + triangle_field tri = triangle_field_polygonise( + density, ctx, 1.0, mtetras); // mtetras or mcubes - char *stl_fname = stl_create_fname(prefix, "iso_surface_tetras", it); + char *stl_fname = + stl_create_fname(prefix, "iso_surface_tetras", 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); - tri = triangle_field_polygonise(density, ctx, 1.0, mcubes); // mtetras or mcubes + tri = triangle_field_polygonise( + density, ctx, 1.0, mcubes); // mtetras or mcubes stl_fname = stl_create_fname(prefix, "iso_surface_cubes", it); - stl = stl_file_create(stl_fname); + stl = stl_file_create(stl_fname); stl_file_write_triangle_field(&stl, tri); stl_file_close(&stl); free(stl_fname); triangle_field_deallocate(tri); - } - f = lbm_collide_and_stream(ctx, &f, vel_ini, d, units.omega, delta_iter); + f = lbm_collide_and_stream( + ctx, &f, vel_ini, d, units.omega, delta_iter); // TODO add version with rho/u as output - // f = lbm_collide_and_stream_rho_u(ctx, &f, &f_rho, &f_vel, units.omega, delta_iter); + // f = lbm_collide_and_stream_rho_u(ctx, &f, &f_rho, &f_vel, + // units.omega, delta_iter); i += 1; } - + futhark_context_sync(ctx); clock_gettime(CLOCK_MONOTONIC, &finish); T elapsed = (finish.tv_sec - start.tv_sec) * 1e6; elapsed += (finish.tv_nsec - start.tv_nsec) / 1e3; printf("elapsed microseconds = %f\n", elapsed); - printf("approximative MLUSP = %f\n", (T)units.nx * (T)units.ny * (T)units.nz * (T)i * (T)delta_iter / elapsed); + printf("approximative MLUSP = %f\n", (T)units.nx * (T)units.ny * + (T)units.nz * (T)i * + (T)delta_iter / elapsed); - if (prefix != NULL) { printf("output at iter %d\n", iter); tensor_field_from_futhark_f32_4d_inplace(f, ctx, lattice); @@ -152,17 +171,18 @@ int main(int argc, char **argv) { char *fname = vtk_create_fname(prefix, "out", iter); - box domain = box_create(0, lattice.nx-1, 0, lattice.ny-1, 0, lattice.nz-1); + box domain = + box_create(0, lattice.nx - 1, 0, lattice.ny - 1, 0, lattice.nz - 1); vtk_file vtk = vtk_file_create_default(fname, domain, 1.0); vtk_file_write_scalar_field(&vtk, &density, domain, "density", base64); - vtk_file_write_tensor_field(&vtk, &velocity, domain, "velocity", base64); + vtk_file_write_tensor_field( + &vtk, &velocity, domain, "velocity", base64); vtk_file_close(&vtk); - + free(fname); } - tensor_field_deallocate(lattice); tensor_field_deallocate(velocity); scalar_field_deallocate(density); @@ -171,12 +191,12 @@ int main(int argc, char **argv) { futhark_free_4d(ctx, f); futhark_free_4d(ctx, d); futhark_free_4d(ctx, vel_ini); -// TODO add version with rho/u output of futhark + // TODO add version with rho/u output of futhark // futhark_free_4d(ctx, f_vel); // futhark_free_f32_3d(ctx, f_rho); - futhark_context_config_free(cfg); futhark_context_free(ctx); + futhark_context_config_free(cfg); return EXIT_SUCCESS; }