ArrayFire: pde/boltzmann_cfd.cpp (original) (raw)
#include
#include
#include
static const float ex_vals[] = {1.0, 1.0, 1.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0};
static const float ey_vals[] = {1.0, 0.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, -1.0};
static const float wt_vals[] = {1.0f / 36.0f, 4.0f / 36.0f, 1.0f / 36.0f,
4.0f / 36.0f, 16.0f / 36.0f, 4.0f / 36.0f,
1.0f / 36.0f, 4.0f / 36.0f, 1.0f / 36.0f};
static const int opposite_indices[] = {8, 7, 6, 5, 4, 3, 2, 1, 0};
struct Simulation {
size_t grid_width;
size_t grid_height;
float density;
float velocity;
float reynolds;
};
Simulation create_simulation(uint32_t grid_width, uint32_t grid_height,
float density, float velocity, float reynolds,
const char* ux_image_filename,
const char* uy_image_filename,
const char* boundaries_filename) {
Simulation sim;
sim.grid_width = grid_width;
sim.grid_height = grid_height;
sim.velocity = velocity;
sim.density = density;
sim.reynolds = reynolds;
try {
std::cerr << e.what() << std::endl;
sim.ux = af::constant(0, grid_width, grid_height, 3);
}
auto ux_dim = sim.ux.dims();
if (ux_dim[0] != grid_width || ux_dim[1] != grid_height) {
std::cerr
<< "Fluid flow ux image has dimensions different to the simulation"
<< std::endl;
throw std::runtime_error{
"Fluid flow ux image has dimensions different to the simulation"};
}
try {
std::cerr << e.what() << std::endl;
sim.uy = af::constant(0, grid_width, grid_height, 3);
}
auto uy_dim = sim.uy.dims();
if (uy_dim[0] != grid_width || uy_dim[1] != grid_height) {
std::cerr
<< "Fluid flow uy image has dimensions different to the simulation"
<< std::endl;
throw std::runtime_error{
"Fluid flow uy image has dimensions different to the simulation"};
}
try {
sim.set_boundaries = af::loadImage(boundaries_filename, false);
std::cerr << e.what() << std::endl;
sim.set_boundaries = af::constant(0, grid_width, grid_height);
}
auto b_dim = sim.set_boundaries.dims();
if (b_dim[0] != grid_width || b_dim[1] != grid_height) {
std::cerr
<< "Fluid boundary image has dimensions different to the simulation"
<< std::endl;
throw std::runtime_error{
"Fluid boundary image has dimensions different to the simulation"};
}
velocity / 255.f;
velocity / 255.f;
sim.set_boundaries = sim.set_boundaries.T() > 0;
return sim;
}
void initialize(Simulation& sim) {
auto& ux = sim.ux;
auto& uy = sim.uy;
auto& rho = sim.rho;
auto& sigma = sim.sigma;
auto& f = sim.f;
auto& feq = sim.feq;
auto& ex = sim.ex;
auto& ey = sim.ey;
auto& wt = sim.wt;
auto& ex_ = sim.ex_;
auto& ey_ = sim.ey_;
auto& ex_T = sim.ex_T;
auto& ey_T = sim.ey_T;
auto& wt_T = sim.wt_T;
auto density = sim.density;
auto velocity = sim.velocity;
auto xcount = sim.grid_width;
auto ycount = sim.grid_height;
ex_ = af::tile(ex, xcount, ycount, 1);
ey_ = af::tile(ey, xcount, ycount, 1);
auto edotu = ex_ * ux + ey_ * uy;
auto udotu = ux * ux + uy * uy;
feq = rho * wt *
((edotu * edotu * 4.5f) - (udotu * 1.5f) + (edotu * 3.0f) + 1.0f);
f = feq;
}
void collide_stream(Simulation& sim) {
auto& ux = sim.ux;
auto& uy = sim.uy;
auto& rho = sim.rho;
auto& sigma = sim.sigma;
auto& f = sim.f;
auto& feq = sim.feq;
auto& set_boundaries = sim.set_boundaries;
auto& ex = sim.ex;
auto& ey = sim.ey;
auto& wt = sim.wt;
auto& ex_ = sim.ex_;
auto& ey_ = sim.ey_;
auto& ex_T = sim.ex_T;
auto& ey_T = sim.ey_T;
auto& wt_T = sim.wt_T;
auto density = sim.density;
auto velocity = sim.velocity;
auto reynolds = sim.reynolds;
auto xcount = sim.grid_width;
auto ycount = sim.grid_height;
const float viscosity =
velocity * std::sqrt(static_cast<float>(xcount * ycount)) / reynolds;
const float tau = 0.5f + 3.0f * viscosity;
const float csky = 0.16f;
auto edotu = ex_ * ux + ey_ * uy;
auto udotu = ux * ux + uy * uy;
feq =
rho * wt * (edotu * edotu * 4.5f - udotu * 1.5f + edotu * 3.0f + 1.0f);
auto taut =
af::sqrt(sigma * (csky * csky * 18.0f * 0.25f) + (tau * tau * 0.25f)) -
(tau * 0.5f);
auto fplus = f - (f - feq) / (taut + tau);
for (int i = 0; i < 9; ++i) {
int xshift = static_cast<int>(ex_vals[i]);
int yshift = static_cast<int>(ey_vals[i]);
}
f = fplus;
ux_top =
ux_bot =
uy_top =
uy_bot =
auto ux_rht = af::tile(ux.cols(ycount - 3, ycount - 1), af::dim4(1, 3));
auto uy_rht = af::tile(uy.cols(ycount - 3, ycount - 1), af::dim4(1, 3));
auto ubdoute_top = ux_top * ex_T + uy_top * ey_T;
auto ubdoute_bot = ux_bot * ex_T + uy_bot * ey_T;
auto ubdoute_lft = ux_lft * ex_T + uy_lft * ey_T;
auto ubdoute_rht = ux_rht * ex_T + uy_rht * ey_T;
6.0 * density * wt_T * ubdoute_top;
6.0 * density * wt_T * ubdoute_bot;
6.0 * density * wt_T * ubdoute_lft;
6.0 * density * wt_T * ubdoute_rht;
for (int i = 0; i < 9; ++i) {
int xshift = static_cast<int>(ex_vals[i]);
int yshift = static_cast<int>(ey_vals[i]);
if (xshift == 1)
if (xshift == -1)
f(xcount - 2, af::span, opposite_indices[i]) =
if (yshift == 1)
if (yshift == -1)
f(af::span, ycount - 2, opposite_indices[i]) =
}
}
void update(Simulation& sim) {
auto& ux = sim.ux;
auto& uy = sim.uy;
auto& rho = sim.rho;
auto& sigma = sim.sigma;
auto& f = sim.f;
auto& feq = sim.feq;
auto& ex = sim.ex;
auto& ey = sim.ey;
auto result = af::sum(f * e_tile, 2);
result /= rho;
auto e_product = af::join(3, ex * ex, ex * ey * std::sqrt(2), ey * ey);
}
af::array generate_image(size_t width, size_t height, const Simulation& sim) {
const auto& ux = sim.ux;
const auto& uy = sim.uy;
const auto& boundaries = sim.set_boundaries;
auto velocity = sim.velocity;
float image_scale =
static_cast<float>(width) / static_cast<float>(sim.grid_width - 1);
auto val = af::sqrt(ux * ux + uy * uy) / velocity;
af::replace(val, val != 0 || !boundaries, -1.0);
if (width != sim.grid_width || height != sim.grid_height)
val =
val = val.T();
auto image2 = image;
auto tile_val = af::tile(val, 1, 1, 3);
return image;
}
void lattice_boltzmann_cfd_demo() {
const size_t len = 128;
const size_t grid_width = len;
const size_t grid_height = len;
int height = static_cast<int>(grid_width * scale);
int width = static_cast<int>(grid_height * scale);
af::Window window(height, width, "Driven Cavity Flow");
int frame_count = 0;
int max_frames = 20000;
int simulation_frames = 100;
float total_time = 0;
float total_time2 = 0;
const float density = 2.7f;
const float velocity = 0.35f;
const float reynolds = 1e5f;
const char* ux_image = ASSETS_DIR "/examples/images/default_ux.bmp";
const char* uy_image = ASSETS_DIR "/examples/images/default_uy.bmp";
const char* set_boundary_image =
ASSETS_DIR "/examples/images/default_boundary.bmp";
{
}
{
}
Simulation sim =
create_simulation(grid_width, grid_height, density, velocity, reynolds,
ux_image, uy_image, set_boundary_image);
initialize(sim);
while (!window.close() && frame_count != max_frames) {
auto begin = std::chrono::high_resolution_clock::now();
collide_stream(sim);
update(sim);
auto end = std::chrono::high_resolution_clock::now();
auto duration =
std::chrono::duration_caststd::chrono::microseconds(end - begin)
.count();
total_time += duration;
total_time2 += duration * duration;
if (frame_count % simulation_frames == 0) {
auto image = generate_image(width, height, sim);
window.image(image);
float avg_time = total_time / (float)simulation_frames;
float stdv_time = std::sqrt(total_time2 * simulation_frames -
total_time * total_time) /
(float)simulation_frames;
std::cout << "Average Simulation Step Time: (" << avg_time
<< " +/- " << stdv_time
<< ") us; Total simulation time: " << total_time
<< " us; Simulation Frames: " << simulation_frames
<< std::endl;
total_time = 0;
total_time2 = 0;
}
frame_count++;
}
}
int main(int argc, char** argv) {
int device = argc > 1 ? std::atoi(argv[1]) : 0;
try {
std::cout << "** ArrayFire CFD Simulation Demo\n\n";
lattice_boltzmann_cfd_demo();
std::cerr << e.what() << std::endl;
return -1;
}
return 0;
}
Window object to render af::arrays.
A multi dimensional data container.
Generic object that represents size and shape.
An ArrayFire exception class.
virtual const char * what() const
Returns an error message for the exception in a string format.
@ f32
32-bit floating point values
AFAPI array pow(const array &base, const array &exponent)
C++ Interface to raise a base to a power (or exponent).
AFAPI array sqrt(const array &in)
C++ Interface to evaluate the square root.
array::array_proxy rows(int first, int last)
Returns a reference to sequence of rows.
array constant(T val, const dim4 &dims, const dtype ty=(af_dtype) dtype_traits< T >::ctype)
C++ Interface to generate an array with elements set to a specified value.
AFAPI array iota(const dim4 &dims, const dim4 &tile_dims=dim4(1), const dtype ty=f32)
C++ Interface to generate an array with [0, n-1] values modified to specified dimensions and tiling.
AFAPI void replace(array &a, const array &cond, const array &b)
C++ Interface to replace elements of an array with elements of another array.
AFAPI void setDevice(const int device)
Sets the current device.
AFAPI void sync(const int device=-1)
Blocks until the device is finished processing.
AFAPI array loadImage(const char *filename, const bool is_color=false)
C++ Interface for loading an image.
AFAPI array join(const int dim, const array &first, const array &second)
C++ Interface to join 2 arrays along a dimension.
AFAPI array moddims(const array &in, const dim4 &dims)
C++ Interface to modify the dimensions of an input array to a specified shape.
AFAPI array shift(const array &in, const int x, const int y=0, const int z=0, const int w=0)
C++ Interface to shift an array.
AFAPI array tile(const array &in, const unsigned x, const unsigned y=1, const unsigned z=1, const unsigned w=1)
C++ Interface to generate a tiled array.
AFAPI array product(const array &in, const int dim=-1)
C++ Interface to multiply array elements over a given dimension.
AFAPI array sum(const array &in, const int dim=-1)
C++ Interface to sum array elements over a given dimension.
AFAPI array approx2(const array &in, const array &pos0, const array &pos1, const interpType method=AF_INTERP_LINEAR, const float off_grid=0.0f)
C++ Interface for data interpolation on two-dimensional signals.
AFAPI array scale(const array &in, const float scale0, const float scale1, const dim_t odim0=0, const dim_t odim1=0, const interpType method=AF_INTERP_NEAREST)
C++ Interface for scaling an image.
AFAPI int end
A special value representing the last value of an axis.
AFAPI seq span
A special value representing the entire axis of an af::array.