# FluidX3D Documentation - How to get started?
## 0. Install GPU Drivers and OpenCL Runtime
(click to expand section)
- **Windows**
GPUs
- Download and install the [AMD](https://www.amd.com/en/support/download/drivers.html)/[Intel](https://www.intel.com/content/www/us/en/download/785597/intel-arc-iris-xe-graphics-windows.html)/[Nvidia](https://www.nvidia.com/Download/index.aspx) GPU Drivers, which contain the OpenCL Runtime.
- Reboot.
CPUs
- Download and install the [Intel CPU Runtime for OpenCL](https://www.intel.com/content/www/us/en/developer/articles/technical/intel-cpu-runtime-for-opencl-applications-with-sycl-support.html) (works for both AMD/Intel CPUs).
- Reboot.
AMD GPUs
- Download and install [AMD GPU Drivers](https://www.amd.com/en/support/download/linux-drivers.html), which contain the OpenCL Runtime, with:
```bash
sudo apt update && sudo apt upgrade -y
sudo apt install -y g++ git make ocl-icd-libopencl1 ocl-icd-opencl-dev
mkdir -p ~/amdgpu
wget -P ~/amdgpu https://repo.radeon.com/amdgpu-install/25.35.1/ubuntu/noble/amdgpu-install_7.2.1.70201-1_all.deb
sudo apt install -y ~/amdgpu/amdgpu-install*.deb
sudo amdgpu-install -y --usecase=graphics,rocm,opencl --opencl=rocr
sudo usermod -a -G render,video $(whoami)
rm -r ~/amdgpu
sudo shutdown -r now
```
Intel GPUs
- Intel GPU Drivers come already installed since Linux Kernel 6.2, but they don't contain the OpenCL Runtime.
- The the [OpenCL Runtime](https://github.com/intel/compute-runtime/releases) has to be installed separately with:
```bash
sudo apt update && sudo apt upgrade -y
sudo apt install -y g++ git make ocl-icd-libopencl1 ocl-icd-opencl-dev intel-opencl-icd
sudo usermod -a -G render $(whoami)
sudo shutdown -r now
```
Nvidia GPUs
- Download and install [Nvidia GPU Drivers](https://www.nvidia.com/Download/index.aspx), which contain the OpenCL Runtime, with:
```bash
sudo apt update && sudo apt upgrade -y
sudo apt install -y g++ git make ocl-icd-libopencl1 ocl-icd-opencl-dev nvidia-driver-580
sudo shutdown -r now
```
CPUs
- Option 1: Download and install the [oneAPI DPC++ Compiler](https://github.com/intel/llvm/releases?q=%22oneAPI+DPC%2B%2B+Compiler+dependencies%22) and [oneTBB](https://github.com/uxlfoundation/oneTBB/releases) with:
```bash
export OCLV="oclcpuexp-2025.21.10.0.10_160000_rel"
export TBBV="oneapi-tbb-2023.0.0"
sudo apt update && sudo apt upgrade -y
sudo apt install -y g++ git make ocl-icd-libopencl1 ocl-icd-opencl-dev
sudo mkdir -p ~/cpurt /opt/intel/${OCLV} /etc/OpenCL/vendors /etc/ld.so.conf.d
sudo wget -P ~/cpurt https://github.com/intel/llvm/releases/download/2025-WW45/${OCLV}.tar.gz
sudo wget -P ~/cpurt https://github.com/uxlfoundation/oneTBB/releases/download/v2023.0.0/${TBBV}-lin.tgz
sudo tar -zxvf ~/cpurt/${OCLV}.tar.gz -C /opt/intel/${OCLV}
sudo tar -zxvf ~/cpurt/${TBBV}-lin.tgz -C /opt/intel
echo /opt/intel/${OCLV}/x64/libintelocl.so | sudo tee /etc/OpenCL/vendors/intel_expcpu.icd
echo /opt/intel/${OCLV}/x64 | sudo tee /etc/ld.so.conf.d/libintelopenclexp.conf
sudo ln -sf /opt/intel/${TBBV}/lib/intel64/gcc4.8/libtbb.so /opt/intel/${OCLV}/x64
sudo ln -sf /opt/intel/${TBBV}/lib/intel64/gcc4.8/libtbbmalloc.so /opt/intel/${OCLV}/x64
sudo ln -sf /opt/intel/${TBBV}/lib/intel64/gcc4.8/libtbb.so.12 /opt/intel/${OCLV}/x64
sudo ln -sf /opt/intel/${TBBV}/lib/intel64/gcc4.8/libtbbmalloc.so.2 /opt/intel/${OCLV}/x64
sudo ldconfig -f /etc/ld.so.conf.d/libintelopenclexp.conf
sudo rm -r ~/cpurt
```
- Option 2: Download and install [PoCL](https://portablecl.org/) with:
```bash
sudo apt update && sudo apt upgrade -y
sudo apt install -y g++ git make ocl-icd-libopencl1 ocl-icd-opencl-dev pocl-opencl-icd
```
ARM GPUs
- Download the [Termux `.apk`](https://github.com/termux/termux-app/releases) and install it.
- In the Termux app, run:
```bash
apt update && apt upgrade -y
apt install -y clang git make
```
## 1. Download FluidX3D
- [Download](https://github.com/ProjectPhysX/FluidX3D/archive/refs/heads/master.zip) and unzip the source code, or clone with:
```bash
git clone https://github.com/ProjectPhysX/FluidX3D.git && cd FluidX3D
```
- To update FluidX3D:
- Make a backup of your changes.
- Run:
```bash
git stash
git pull origin master
git stash pop
```
## 2. Compiling the Source Code
- There is no "installation" of FluidX3D. Instead, you have to compile the source code yourself.
- I have made this as easy as possible and this documentation will guide you through it. Nonetheless, some basic programming experience with C++ would be good for the setup scripts.
- First, compile the code as-is; this is the standard FP32 benchmark test case. By default, the fastest installed GPU will be selected automatically. Compile time is about 5 seconds.
### Windows
- Download and install [Visual Studio Community](https://visualstudio.microsoft.com/de/vs/community/). In Visual Studio Installer, add:
- Desktop development with C++
- MSVC v142
- Windows 10 SDK
- Open [`FluidX3D.sln`](FluidX3D.sln) in [Visual Studio Community](https://visualstudio.microsoft.com/de/vs/community/).
- Compile and run by clicking the ► Local Windows Debugger button.
- To select a specific GPU, open Windows CMD in the `FluidX3D` folder (type `cmd` in File Explorer in the directory field and press Enter), then run `bin\FluidX3D.exe 0` to select device `0`. You can also select multiple GPUs with `bin\FluidX3D.exe 0 1 3 6` if the setup is [configured as multi-GPU](#the-lbm-class).
### Linux / macOS / Android
- Compile and run with:
```bash
chmod +x make.sh
./make.sh
```
- Compiling requires [`g++`](https://gcc.gnu.org/) with `C++17`, which is supported since version `8` (check with `g++ --version`). If you have [`make`](https://www.gnu.org/software/make/) installed (check with `make --version`), compiling will will be faster using multiple CPU cores; otherwise compiling falls back to using a single CPU core.
- To select a specific GPU, enter `./make.sh 0` to compile+run, or `bin/FluidX3D 0` to run on device `0`. You can also select multiple GPUs with `bin/FluidX3D 0 1 3 6` if the setup is [configured as multi-GPU](#the-lbm-class).
- Operating system (Linux/macOS/Android) and X11 support (required for [`INTERACTIVE_GRAPHICS`](src/defines.hpp)) are detected automatically. In case problems arise, you can still manually select [`target=...`](make.sh#L13) in [`make.sh`](make.sh#L13).
- On macOS, [`INTERACTIVE_GRAPHICS`](src/defines.hpp) mode requires [XQuartz](https://www.xquartz.org/) to be installed.
- On Android, [`INTERACTIVE_GRAPHICS`](src/defines.hpp) mode is not supported, as no X11 is available. You can still use [`INTERACTIVE_GRAPHICS_ASCII`](src/defines.hpp) though, or [render video](#video-rendering) to the hard drive with regular [`GRAPHICS`](src/defines.hpp) mode.
## 3. Go through Sample Setups
- Now open [`src/setup.cpp`](src/setup.cpp). In here are all the sample setups, each one being a `void main_setup() {...}` function block written in C++. Uncomment one of them, maybe start top-to-bottom.
- In the line where the `main_setup()` function starts, it says "required extensions in defines.hpp:", followed by a list of extensions in capital letters. Head over to [`src/defines.hpp`](src/defines.hpp) and comment out
```c
//#define BENCHMARK
```
with a `//`. Then, uncomment all of the extensions required for the setup by removing the `//` in front of the corresponding line.
- Finally, [compile](#2-compiling-the-source-code) and run the setup with the ► Local Windows Debugger button (Windows) or `./make.sh` (Linux/macOS/Android).
- Once the interactive graphics window opens, press key P to start/pause the simulation, and press H to show the help menu for keyboard controls and visualization settings.
- Go through some of the sample setups this way, get familiar with their code structure and test the graphics mode.
## 4. Keyboard/Mouse Controls for [`INTERACTIVE_GRAPHICS`](src/defines.hpp)/[`_ASCII`](src/defines.hpp)
| Key | Function |
| :-----------------------: | :--------------------------------------------------------------------------------------------------------------------------------- |
| | |
| P | start/pause the simulation |
| H | show/hide help menu for keyboard controls and visualization settings |
| Esc
Alt+F4 | quit |
| | |
| Mouse
I
J K L | rotate camera |
| Scrollwheel
+ - | zoom (centered camera mode) or camera movement speed (free camera mode) |
| Mouseclick
U | toggle rotation with Mouse and angle snap rotation with I J K L |
| F | toggle centered/free camera mode |
| W
A S D
Space C | move free camera |
| Y X | adjust camera field of view |
| R | toggle camera autorotation |
| G | print current camera position/rotation in console as copy/paste command |
| V | toggle stereoscopic rendering for VR |
| B | toggle VR-goggles/3D-TV mode for stereoscopic rendering |
| N M | adjust eye distance for stereoscopic rendering |
| | |
| 1 | flag wireframe / solid surface (and force vectors on solid cells or surface pressure if the extension is used) |
| 2 | velocity field |
| 3 | streamlines |
| 4 | vorticity (velocity-colored Q-criterion isosurface) |
| 5 | rasterized free surface |
| 6 | raytraced free surface |
| 7 | particles |
| T | toggle slice visualization mode |
| Z | toggle field visualization mode |
| Q E | move slice in slice visualization mode |
## 5. Writing your own Setups
### The LBM Class
- For initializing the simulation box, use call
```c
LBM lbm(Nx, Ny, Nz, nu, ...);
```
constructor. `Nx`/`Ny`/`Nz` is the grid resolution and `nu` is the kinematic shear viscosity in [LBM units](#unit-conversion).
- To use multiple GPUs, use
```c
LBM lbm(Nx, Ny, Nz, Dx, Dy, Dz, nu, ...);
```
with `Dx`/`Dy`/`Dz` indicating how many domains (GPUs) there are in each spatial direction. The product `Dx`×`Dy`×`Dz` is the total number of domains (GPUs).
- As long as the `lbm` object is in scope, you can access the memory. As soon as it goes out of scope, all memory associated with the current simulation is freed again.
- The grid resolution `Nx`/`Ny`/`Nz` ultimately determines the VRAM occupation. Quite often it's not obvious at which resolution you'll overshoot the VRAM capacity of the GPU(s). To aid with this, there is the function:
```c
const uint3 lbm_N = resolution(float3(1.0f, 2.0f, 0.5f), 2000u);
```
This takes as inputs the desired aspect ratio of the simulation box and the VRAM occupation in MB, and returns the grid resolution as a `uint3` with `.x`/`.y`/`.z` components. You can also directly feed the `uint3` into the LBM constructor as resolution:
```c
LBM lbm(lbm_N, nu, ...);
```
### Unit Conversion
- The LBM simulation uses a different unit system from SI units, where density `rho=1` and velocity `u≈0.001-0.1`, because floating-point arithmetic is most accurate close to `1`.
- To ease unit conversion from SI to LBM units and back, there is the [`units.hpp`](src/units.hpp) struct. By calling
```c
units.set_m_kg_s(lbm_length, lbm_velocity, lbm_density=1, si_length, si_velocity, si_density);
```
the base unit conversion factors [m], [kg], [s] are calculated and stored in the `units` struct. Thereafter, any of the conversion functions from [`src/units.hpp`](src/units.hpp) can be used to go from SI to LBM units and back, such as `lbm_nu = units.nu(si_nu)` to convert the kinematic viscosity from SI to LBM units.
- A good beginner example for this is the "[aerodynamics of a cow](src/setup.cpp)" setup.
### Initial and Boundary Conditions
- If not explicitly set, by default all cells have the default values `rho=1`, `u=0`, `flags=0`.
- The initial/boundary conditions of single grid cells are set in a parallelized loop that iterates over the entire grid:
```c
const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); parallel_for(lbm.get_N(), [&](ulong n) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
// ...
});
```
Within this loop, you can set the density, velocity and flags of each cell individually by assigning values to `lbm.rho[n]`, `lbm.u.x[n]`, `lbm.u.y[n]`, `lbm.u.z[n]` and `lbm.flags[n]`. The `n` here is the linearized 3D grid index, corresponding to an (`x`|`y`|`z`) position via the function `lbm.coordinates(n, x, y, z)`.
- For example, to set solid boundaries on all sides of the simulation box, within the initialization loop, write:
```c
if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) {
lbm.flags[n] = TYPE_S;
}
```
- For primitive geometry, such as spheres, ellipsoids, cubes, cuboids, cylinders, codes, pipes, triangles, inclined planes, or toruses can be set with the functions from [`shapes.hpp`](src/shapes.hpp). Example to insert a cylinder:
```c
if(cylinder(x, y, z, lbm.center(), float3(axis_x, axis_z, axis_z), radius) {
lbm.flags[n] = TYPE_S;
}
```
- Available Boundary Conditions
- Periodic Boundaries
- All box sides where no solid (`TYPE_S`) or other boundary type are set will remain periodic boundaries.
- If strict mass conservation is required (for example flow through a linear pipe), use periodic boundaries (i.e. don't set any boundary type on the cells at these simulation box sides), and drive the flow with a volume force (equivalent to a pressure gradient). For this you need to enable (uncomment) the [`VOLUME_FORCE`](src/defines.hpp) extension, and in the [LBM constructor](#the-lbm-class) set the force per volume (`fx`|`fy`|`fz`):
```c
LBM lbm(Nx, Ny, Nz, nu, fx, fy, fz);
```
These force per volume values should not exceed `0.001` in magnitude.
- Solid Boundaries (`TYPE_S`)
- Always available without further extensions.
- Solid boundaries function as solid non-moving objects, and are implemented as non-moving no-slip mid-grid bounce-back boundaries.
- "non-moving" - velocity of the boundary is 0
- "no-slip bounce-back" - flow velocity directly at the boundary is 0 (no-slip condition)
- "mid-grid" - boundary is located exactly in the middle between the boundary cells and adjacent fluid cells
- Usage:
```c
if() {
lbm.flags[n] = TYPE_S;
}
```
- Equilibrium Boundaries (`TYPE_E` with specified density and/or velocity)
- You need to enable (uncomment) the [`EQUILIBRIUM_BOUNDARIES`](src/defines.hpp) extension.
- Equilibrium boundaries function as inflow/outflow boundaries. `lbm.flags[n] = TYPE_E` cells enforce the density `lbm.rho[n]` and/or velocity `lbm.u.x[n]`/`lbm.u.y[n]`/`lbm.u.z[n]` specified during initialization and absorb any incoming shockwaves. If no density/velocity values are specified, they will enforce the default values (`rho=1`, `u=0`).
- Usage: For the specific inflow/outflow cells, set the flag `lbm.flags[n] = TYPE_E` and on the same cells specify either a density `lbm.rho[n]` unequal to `1` or a velocity `lbm.u.x[n]`/`lbm.u.y[n]`/`lbm.u.z[n]` unequal to `0`, or a combination of both.
```c
if() {
lbm.flags[n] = TYPE_E;
lbm.rho[n] = <...>; // TYPE_E cells will retain this fixed density
lbm.u.x[n] = <...>; // TYPE_E cells will retain this fixed velocity
lbm.u.y[n] = <...>; // TYPE_E cells will retain this fixed velocity
lbm.u.z[n] = <...>; // TYPE_E cells will retain this fixed velocity
}
```
- Moving Solid Boundaries (`TYPE_S` with non-zero velocity)
- You need to enable (uncomment) the [`MOVING_BOUNDARIES`](src/defines.hpp) extension.
- Moving solid boundaries can function as solid moving objects, or they can function as inflow/outflow with fixed flow rate (similar to a hydraulic piston). They reflect any incoming shockwaves.
- Moving solid boundaries are implemented as moving no-slip mid-grid bounce-back boundaries.
- "moving" - velocity of the boundary is the specified non-zero value (Dirichlet boundary condition for velocity).
- "no-slip bounce-back" - flow velocity directly at the boundary is equal to the non-zero boundary velocity (no-slip condition)
- "mid-grid" - boundary is located exactly in the middle between the boundary cells and adjacent fluid cells
- Usage: For the specific moving solid cells, set the flag `lbm.flags[n] = TYPE_S` and on the same cells specify a velocity `lbm.u.x[n]`/`lbm.u.y[n]`/`lbm.u.z[n]` unequal to `0`.
```c
if() {
lbm.flags[n] = TYPE_S;
lbm.u.x[n] = <...>; // TYPE_S cells will retain this fixed velocity
lbm.u.y[n] = <...>; // TYPE_S cells will retain this fixed velocity
lbm.u.z[n] = <...>; // TYPE_S cells will retain this fixed velocity
}
```
- You can also use moving solid boundaries as inflow with fixed flow rate with the [`SURFACE`](#surface-extension) extension. Here, in addition to specifying the velocity in `TYPE_S` cells, you need to initialize the layer of cells next to the `TYPE_S` cells with `TYPE_F` flag.
### Running the Simulation
- Call `lbm.run()` (without input parameter, it's infinite time steps) to initialize and execute the setup, or `lbm.run(time_steps)` to execute only a specific number of time steps.
- If you have a [more complicated simulation loop](#video-rendering) where you periodically compute time steps and render images for a video or export data, don't forget to place an `lbm.run(0u)` before that loop. This copies the initial/boundary conditions from CPU RAM to GPU VRAM and initializes the simulation on the GPU, without computing a time step. Without initialization, there is no data in VRAM yet for rendering.
### Loading .stl Files
- For more complex geometries, you can load `.stl` triangle meshes and voxelize them to the Cartesian simulation grid on the GPU(s).
- Create a `FluidX3D/stl/` folder next to the `FluidX3D/src/` folder and download the geometry from websites like [Thingiverse](https://www.thingiverse.com/), or create your own.
- Only binary `.stl` files are supported. Meshes must be watertight (no holes) and all triangles must be oriented such that their normals point to the outside. For conversion from other formats or for splitting composite geometries like helicopter hull and rotors, I recommend [Microsoft 3D Builder](https://apps.microsoft.com/detail/9wzdncrfj3t6) (to download: paste `apps.microsoft.com/detail/9wzdncrfj3t6` in [here](https://store.rg-adguard.net/) and click on `Microsoft.3DBuilder_20.0.4.0_neutral_~_8wekyb3d8bbwe.appxbundle`) on Windows or [Blender](https://www.blender.org/) on Windows/Linux.
- Load and voxelize simple `.stl` files directly with
```c
lbm.voxelize_stl(get_exe_path()+"../stl/mesh.stl", center, rotation, size);
```
This automatically repositions/rescales the mesh to the specified center. Use `lbm.center()` for the simulation box center, or add an offset with a `+float3(offset_x, offset_y, offset_z)`. You can generate and multiply together a rotation matrix like this (example: rotation around the z-axis by 180°, then around the x-axis by 90°):
```c
float3x3 rotation = float3x3(float3(1, 0, 0), radians(90.0f))*float3x3(float3(0, 0, 1), radians(180.0f));
```
- To load composite geometries with several parts without automatic mesh repositioning/rescaling, use
```c
Mesh* mesh_1 = read_stl(const string& path, const float scale=1.0f, const float3x3& rotation=float3x3(1.0f), const float3& offset=float3(0.0f)); // load mesh without automatic repositioning/rescaling
Mesh* mesh_2 = read_stl(const string& path, const float scale=1.0f, const float3x3& rotation=float3x3(1.0f), const float3& offset=float3(0.0f));
mesh_1->scale(const float scale); // manually scale meshes
mesh_2->scale(const float scale);
mesh_1->translate(const float3& translation); // manually reposition meshes
mesh_2->translate(const float3& translation);
lbm.voxelize_mesh_on_device(mesh_1); // voxelize meshes on GPU
lbm.voxelize_mesh_on_device(mesh_2);
```
to load the meshes from the `.stl` files, manually scale/reposition all parts of the mesh the same time, and finally voxelize them on the GPU.
- To aid with repositioning the mesh, there is `lbm.center()` for the center of the simulation box, as well as the min/max bounding-box coordinates of the mesh `mesh->pmin`/`mesh->pmax`, each a `float3` with (`x`|`y`|`z`) components.
- Rotating geometries have to be periodically revoxelized, about every 1-10 LBM time steps. In the main simulation loop in the [`main_setup()`](src/setup.cpp) function, first rotate the triangle mesh, then revoxelize on GPU, then compute a few LBM time steps:
```c
const uint lbm_T = 100000u; // number of LBM time steps to simulate
const uint lbm_dt = 4u; // number of LBM time steps between each mesh revoxelization
lbm.run(0u); // initialize simulation
mesh->set_center(mesh->get_center_of_mass()); // set rotation center of mesh to its center of mass
while(lbm.get_t()rotate(float3x3(float3(0, 0, 1), lbm_omega*(float)lbm_dt)); // rotate the triangle mesh
lbm.voxelize_mesh_on_device(mesh, TYPE_S, center, float3(0.0f), float3(0.0f, 0.0f, lbm_omega)); // revoxelize the rotated triangle mesh, provide the instantaneous angular velocity vector for moving boundaries
lbm.run(lbm_dt); // run lbm_dt LBM time steps
}
```
Here `lbm_omega` is the angular velocity in radians per time step, `lbm_dt` is the number of simulated time steps between revoxelizations, and `float3(0.0f, 0.0f, lbm_omega)` is the instantaneous angular velocity as a vector along the axis of rotation. The largest displacement of the outermost cells should not exceed `1` cell between revoxelizations; set `lbm_omega = lbm_u/lbm_radius` accordingly.
- Have a look at the "[Cessna 172](src/setup.cpp)" and "[Bell 222](src/setup.cpp)" setups for some examples.
### Video Rendering
- For video rendering, disable (comment out) [`INTERACTIVE_GRAPHICS`](src/defines.hpp) and [`INTERACTIVE_GRAPHICS_ASCII`](src/defines.hpp) and enable (uncomment) [`GRAPHICS`](src/defines.hpp) in [`src/defines.hpp`](src/defines.hpp).
- Set the video resolution as [`GRAPHICS_FRAME_WIDTH`](src/defines.hpp)/[`GRAPHICS_FRAME_HEIGHT`](src/defines.hpp) and the background color as [`GRAPHICS_BACKGROUND_COLOR`](src/defines.hpp). You can also adjust the other [`GRAPHICS_...`](src/defines.hpp) options there, such as semi-transparent rendering mode, or adjust the color scale for velocity with [`GRAPHICS_U_MAX`](src/defines.hpp).
- A basic loop for rendering video in the [`main_setup()`](src/setup.cpp) function looks like this:
```c
lbm.graphics.visualization_modes = VIS_FLAG_LATTICE|VIS_Q_CRITERION; // set visualization modes, see all available visualization mode macros (VIZ_...) in defines.hpp
const uint lbm_T = 10000u; // number of LBM time steps to simulate
lbm.run(0u, lbm_T); // initialize simulation
while(lbm.get_t()Mouse to disable mouse rotation, and press G to print the current camera settings as a copy-paste command in the console. Alt+Tab to the console and copy the camera placement command by selecting it with the mouse and right-clicking, then paste it into the [`main_setup()`](src/setup.cpp) function.
- To fly the camera along a smooth path through a list of provided keyframe camera placements, use `catmull_rom` splines:
```c
while(lbm.get_t()<=lbm_T) { // main simulation loop
if(lbm.graphics.next_frame(lbm_T, 30.0f)) {
const float t = (float)lbm.get_t()/(float)lbm_T;
vector camera_positions = {
float3(-0.282220f*(float)Nx, 0.529221f*(float)Ny, 0.304399f*(float)Nz),
float3( 0.806921f*(float)Nx, 0.239912f*(float)Ny, 0.436880f*(float)Nz),
float3( 1.129724f*(float)Nx, -0.130721f*(float)Ny, 0.352759f*(float)Nz),
float3( 0.595601f*(float)Nx, -0.504690f*(float)Ny, 0.203096f*(float)Nz),
float3(-0.056776f*(float)Nx, -0.591919f*(float)Ny, -0.416467f*(float)Nz)
};
vector camera_rx = {
116.0f,
25.4f,
-10.6f,
-45.6f,
-94.6f
};
vector camera_ry = {
26.0f,
33.3f,
20.3f,
25.3f,
-16.7f
};
const float camera_fov = 90.0f;
lbm.graphics.set_camera_free(catmull_rom(camera_positions, t), catmull_rom(camera_rx, t), catmull_rom(camera_ry, t), camera_fov);
lbm.graphics.write_frame(get_exe_path()+"export/");
}
lbm.run(1u, lbm_T);
}
```
- The visualization mode(s) can be specified as `lbm.graphics.visualization_modes` with the [`VIS_...`](src/defines.hpp) macros. You can also set the `lbm.graphics.slice_mode` (`0`=no slice, `1`=x, `2`=y, `3`=z, `4`=xz, `5`=xyz, `6`=yz, `7`=xy) and reposition the slices with `lbm.graphics.slice_x`/`lbm.graphics.slice_y`/`lbm.graphics.slice_z`.
- Exported frames will automatically be assigned the current simulation time step in their name, in the format `bin/export/image-123456789.png`.
- To convert the rendered `.png` images to video, use [FFmpeg](https://ffmpeg.org/):
```bash
ffmpeg -framerate 60 -pattern_type glob -i "export/*/image-*.png" -c:v libx264 -pix_fmt yuv420p -b:v 24M "video.mp4"
```
### Data Export
- At any point in time, you can export volumetric data as binary `.vtk` files with:
```c
lbm.rho.write_device_to_vtk(); // density
lbm.u.write_device_to_vtk(); // velocity
lbm.flags.write_device_to_vtk(); // flags
lbm.F.write_device_to_vtk(); // force, only for FORCE_FIELD extension
lbm.phi.write_device_to_vtk(); // fill fraction, only for SURFACE extension
lbm.T.write_device_to_vtk(); // temperature, only for TEMPERATURE extension
lbm.write_mesh_to_vtk(const Mesh* mesh); // for exporting triangle meshes
```
- These functions first pull the data from the GPU(s) into CPU RAM, and then write it to the hard drive.
- If [unit conversion](#unit-conversion) with `units.set_m_kg_s(...)` was specified, the data in exported `.vtk` files is automaticlally converted to SI units.
- Exported files will automatically be assigned the current simulation time step in their name, in the format `bin/export/u-123456789.vtk`.
- Be aware that these volumetric files can be gigantic in file size, tens of GigaByte for a single file.
- You can view/evaluate the `.vtk` files for example in [ParaView](https://www.paraview.org/).
- It is recommended to use the C++ functionality in the [`main_setup()`](src/setup.cpp) function directly to extract the data of interest and selectively only write that to the hard drive. Therefore, call `lbm.u.read_from_device()` to copy the data from the GPU(s) to CPU RAM, and then you can access it directly, for example
```c
const float lbm_velocity_x = lbm.u.x[lbm.index(x, y, z)];
```
to get the x-velocity at the position (`x`|`y`|`z`) in [LBM units](#unit-conversion).
- To [convert the velocity from LBM to SI units](#unit-conversion), use
```c
const float si_velocity_x = units.si_u(lbm_velocity_x);
```
after having done [unit conversion](#unit-conversion) with `units.set_m_kg_s(...)`.
- You can also export the `.stl` triangle meshes to binary `.vtk` files with:
```c
lbm.write_mesh_to_vtk(const Mesh* mesh);
```
### Lift/Drag Forces
- Enable (uncomment) the [`FORCE_FIELD`](src/defines.hpp) extension. This extension allows computing boundary forces on every solid cell (`TYPE_S`) individually, as well as placing an individual volume force on every fluid cell (not used here).
- In the [`main_setup()`](src/setup.cpp) function, voxelize the mesh with a unique flag combination, such as `(TYPE_S|TYPE_X)` or `(TYPE_S|TYPE_Y)` or `(TYPE_S|TYPE_X|TYPE_Y)`, to distinguish it from all other `(TYPE_S)` cells that might be needed to define other geometry, and compute its center of mass:
```c
lbm.voxelize_mesh_on_device(mesh, TYPE_S|TYPE_X); // voxelize mesh with unique flag combination
const float3 lbm_com = lbm.object_center_of_mass(TYPE_S|TYPE_X); // object center of mass in LBM unit coordinates
```
- To sum over all the individual boundary cells that belong to the object, in the [`main_setup()`](src/setup.cpp) function's main simulation loop call:
```c
const float3 lbm_force = lbm.object_force(TYPE_S|TYPE_X); // force on object
const float3 lbm_torque = lbm.object_torque(lbm_com, TYPE_S|TYPE_X); // torque on object around lbm_com rotation point
```
These functions sum over all cells marked `(TYPE_S|TYPE_X)` that belong to the object. The summation happens GPU-accelerated in VRAM, and only the result is copied to CPU RAM.
- You may also access the force field on individual grid cells. Note that copying the entire `lbm.F` force field from GPU VRAM to CPU RAM is slow:
```c
lbm.F.read_from_device(); // copy entire force field from GPU VRAM to CPU RAM (slow)
lbm_force_x_n = lbm.F.x[lbm.index(x, y, z)]; // access force at one particular grid cell with integer coordinates x, y, z
```
- Finally, [convert from LBM to SI units](#unit-conversion) with
```c
const float si_force_x = units.si_F(lbm_force.x);
```
after having done [unit conversion](#unit-conversion) with `units.set_m_kg_s(...)`.
- See the "[Ahmed body](src/setup.cpp)" sample setup for an example. Note that in the highly turbulent regime, computed body forces are too large by up to a factor 2, because even large resolution is not enough to fully capture the turbulent boundary layer. A wall function is still needed.
## 6. Further LBM Extensions
By now you're already familiar with the [additional boundary types](#initial-and-boundary-conditions) through extensions [`VOLUME_FORCE`](src/defines.hpp), [`FORCE_FIELD`](src/defines.hpp), [`EQUILIBRIUM_BOUNDARIES`](src/defines.hpp), and [`MOVING_BOUNDARIES`](src/defines.hpp). The remaining available model extensions are briefly outlined here:
### [`SURFACE`](src/defines.hpp) Extension
- To simulate free water surfaces, enable (uncomment) the [`SURFACE`](src/defines.hpp) extension.
- All cells then get 3 additional flags: `TYPE_F` (fluid), `TYPE_I` (interface), and `TYPE_G` (gas). Fluid cells are computed with regular LBM. Interface cells account for the extra surface tension forces, if the surface tension coefficient `sigma` is set greater than `0` in the [LBM constructor](#the-lbm-class); the interface is always 1 cell layer thick. Gas cells are not simulated at all and are essentially treated as vacuum.
- If not set otherwise in the [initial conditions](#initial-and-boundary-conditions), all cells are initialized as `TYPE_G` by default. As initial conditions, set all cells that should be fluid to
```c
lbm.flags[n] = TYPE_F;
```
The interface layer will be automatically initialized during initialization with `lbm.run(0u)`.
- Addidionally to the 3 flags, each cell also gets assigned a fill level `lbm.phi[n]`: `1` for fluid cells (`TYPE_F`), `]0,1[` for interface cells (`TYPE_I`), and `0` for gas cells (`TYPE_G`). You can set this fill level at initialization, additionally to the cell flag. Do not forget to set the cell flag. If `lbm.phi[n]` is not set manually, it will automatically be initialized such that all fluid cells get `phi=1`, all interface cells get `phi=0.5`, and all gas clls get `phi=0` assigned.
- For fluid inflow boundary conditions, use [Moving Solid Boundaries](#initial-and-boundary-conditions), and additionally in the adjacent cells to the `TYPE_S` cells with non-zero velocity, set the `TYPE_F` flag.
- See the "[hydraulic jump](src/setup.cpp)" or the "[raindrop impact](src/setup.cpp)" sample setups.
### [`TEMPERATURE`](src/defines.hpp) Extension
- With the [`TEMPERATURE`](src/defines.hpp) extension, FluidX3D can model thermal convection flows. This extension automatically also enables the [`VOLUME_FORCE`](src/defines.hpp) extension.
- In the [LBM constructor](#the-lbm-class), you then need to set the volume force (`fx`|`fy`|`fz`), the thermal diffusion coefficient `alpha`, and the thermal expansion coefficient `beta`, all [in LBM units](#unit-conversion):
```c
LBM lbm(Nx, Ny, Nz, nu, fx, fy, fz, 0.0f, alpha, beta); // the "0.0f" is for the surface tension coefficient sigma which is not used here and has to remain 0
```
- With the extension, each grid cell gets an additional temperature `lbm.T[n]` ([in LBM units](#unit-conversion)) assigned. The default temperature in LBM units is `1`.
- To set temperature boundary conditions, use the flag `TYPE_T` and for the same cells assign a temperature unequal to `1`:
```c
lbm.flags[n] = TYPE_T; // make the cell n a temperature boundary
lbm.T[n] = 1.2f; // set this temperature boundary hotter than average
```
- See the "[Rayleigh-Benard convection](src/setup.cpp)" and "[thermal convection](src/setup.cpp)" setups for two examples.
### [`SUBGRID`](src/defines.hpp) Extension
- Fluid flow is characterized by the Reynolds numberRe = x·u∕nu
with a characteristic length scale `x`, a characteristic velocity `u` and the kinematic shear viscosity `nu`. Larger length scale, larger velocity or smaller viscosity all mean larger Reynolds number.
- The Reynolds number is a unit-less number. A low value Re < 2300 means laminar flow, a high value Re > 2900 means turbulent flow. In between is a transitional regime.
- For very large Reynolds number Re > 100000, the LBM solver becomes [unstable](#7-suitable-parameters-and-simulation-instability), as tiny, very fast rotating vortices can be present in the flow field, and too fast velocity and shear rate makes the simulation blow up.
- To tackle this problem, there is subgrid models that model vortices smaller than single grid cells. This works by increasing the effective vscosity where the shear rate is large and lots of small eddies are assumed to be present. Coincidentally, locations of high shear rate and low viscosity cause instability, so increasing effective viscosity there keeps the simulation stable.
- The subgrid model in FLuidX3D is the Smagorinsky-Lilly model. You can enable it with the [`SUBGRID`](src/defines.hpp) extension.
- There is no additional performance cost for this extension.
### [`PARTICLES`](src/defines.hpp) Extension
- By default, the LBM is a grid-based simulation, so there are no particles.
- But the [`PARTICLES`](src/defines.hpp) extension allows to add particles to the simulation, either as passive tracers or as 2-way-coupled particles that can do floating/sedimentation.
- For passive tracers, only enable the [`PARTICLES`](src/defines.hpp) extension, and in the [LBM constructor](#the-lbm-class) simply add the particle count:
```c
LBM lbm(Nx, Ny, Nz, nu, 50000u); // this will create 50000 particles
```
- Then, in [initialization](#initial-and-boundary-conditions), make a loop over all particles (outside of the initialization loop that iterates over all grid cells):
```c
uint seed = 42u;
for(ulong n=0ull; nlength(); n++) {
lbm.particles->x[n] = random_symmetric(seed, 0.5f*lbm.size().x); // this will palce the particles randomly anywhere in the simulation box
lbm.particles->y[n] = random_symmetric(seed, 0.5f*lbm.size().y);
lbm.particles->z[n] = random_symmetric(seed, 0.5f*lbm.size().z);
}
```
- Note that the position (`0`|`0`|`0`) for particles corresponds to the simulation box center.
- For 2-way-coupled particles, additionally enable the [`VOLUME_FORCE`](src/defines.hpp) and [`FORCE_FIELD`](src/defines.hpp) extensions, and in the [LBM constructor](#the-lbm-class) add the particle density ([in LBM units](#unit-conversion)) unequal to `1`:
```c
LBM lbm(Nx, Ny, Nz, nu, 50000u, 1.2f); // this will create 50000 particles that are more dense than the fluid and will sink to the bottom
```
## 7. Suitable Parameters and Simulation Instability
- Sometimes in the velocity field or streamlines visualization, you will see fuzzyness, or something that looks like a rapidly growing white crystal, blowing up from a certain point and filling the entire simulation box. This is instability, i.e. when velocities turn `NaN` or `Inf`.
- Often times, the cause of instability is an unfortunate choice of unsuitable parameters:
- too high/low density `rho` (ideally should be very close to `1` at all times)
- too high velocity `u` (must never exceed `0.57` anywhere in the box, ideally should be somewhere around `0.075`, but can be as small as `0.001`)
- too low kinematic shear viscosity `nu` (ideally close to `1/6`, becomes unstable when it's very very close to `0` (then enable the [`SUBGRID`](src/defines.hpp) extension), and should not exceed `3`)
- too high force per volume (`fx`|`fy`|`fz`) (should not exceed `0.001` in magnitude)
- too high surface tension coefficient `sigma` (should not exceed `0.1`)
- The best parametrization for LBM simulations is an art in itself and needs some practice.