Interest in surface export functionality?

[nhemsley]

Hello,

I have recently implemented surface export of fluid surfaces to .stl files, it is in a tree with a bunch of other changes, related to a project I am working on. If there is interest in merging this with FluidX3D proper, I can look into extracting the functionality as a 'clean' branch.

This has helped greatly to view, and replay the output of a simulation, in a format, and size which is not so prohibitive. Output from a rtx 5060 5GB simulation, can be replayed in real-time. This can probably be extended to vorticity etc.

Currently it only works on one GPU, as stitching between GPU's requires extra changes. I would like to solve this problem somehow.

As it stands the algorithm does the following (if i remember correctly):

Count the number of triangles, on first frame, copy to host.
Set the triangle buffer to be 1.x multiple of this.
run the existing marching cubes as a kernel function, storing in the surface vertex buffer, with the triangle count.
copy to host, output to stl.
If the triangle count approaches a certain fraction of the triangle buffer, resize the buffer.

Would you be interested in integrating this into FluidX3D?

If so I can look at extracting this functionality in a form easily merged/rebased.

Thanks, Nick

relevant opencl kernel code:

#ifdef SURFACE_EXPORT
//================================================================================================
//                                  SURFACE EXPORT
//================================================================================================
// Functions for exporting fluid surface mesh data for external visualization/processing

kernel void count_surface_triangles(const global float *phi, global uint *triangle_count
#ifdef SURFACE_EXPORT_IGNORE_BOUNDARIES
    , const global uchar *flags
#endif
) {
  const uxx n = get_global_id(0);
  if(n >= def_N) return;

  // Get 3D coordinates
  const uint3 xyz = coordinates(n);
  const uint x = xyz.x;
  const uint y = xyz.y;
  const uint z = xyz.z;

  // Skip boundary cells
  if(x >= def_Nx-1u || y >= def_Ny-1u || z >= def_Nz-1u) return;

  // Get indices of 8 cube corners
  const uxx x0 = (uxx)x, xp = (uxx)(x+1u);
  const uxx y0 = (uxx)(y*def_Nx), yp = (uxx)((y+1u)*def_Nx);
  const uxx z0 = (uxx)(z)*(uxx)(def_Ny*def_Nx), zp = (uxx)(z+1u)*(uxx)(def_Ny*def_Nx);

#ifdef SURFACE_EXPORT_IGNORE_BOUNDARIES
  // Check if any of the 8 cube corners is a solid boundary cell
  const uxx corners[8] = {
    x0+y0+z0, xp+y0+z0, xp+y0+zp, x0+y0+zp,
    x0+yp+z0, xp+yp+z0, xp+yp+zp, x0+yp+zp
  };
  for(uint i = 0u; i < 8u; i++) {
    if(flags[corners[i]] & TYPE_S) return; // Skip this cube if any corner touches a solid boundary
  }
#endif

  float v[8];
  v[0] = phi[x0+y0+z0]; v[1] = phi[xp+y0+z0];
  v[2] = phi[xp+y0+zp]; v[3] = phi[x0+y0+zp];
  v[4] = phi[x0+yp+z0]; v[5] = phi[xp+yp+z0];
  v[6] = phi[xp+yp+zp]; v[7] = phi[x0+yp+zp];

  // Check if cube intersects surface
  float3 triangles[15];
  const uint tn = marching_cubes(v, 0.5f, triangles);

  if(tn > 0u) {
    atomic_add(triangle_count, tn);
  }
}

kernel void export_surface(const global float *phi, global float *vertices,
                          global uint *triangle_count, const ulong max_triangles
#ifdef SURFACE_EXPORT_IGNORE_BOUNDARIES
    , const global uchar *flags
#endif
) {
  const uxx n = get_global_id(0);
  if(n >= def_N) return;

  // Get 3D coordinates
  const uint3 xyz = coordinates(n);
  const uint x = xyz.x;
  const uint y = xyz.y;
  const uint z = xyz.z;

  // Skip boundary cells
  if(x >= def_Nx-1u || y >= def_Ny-1u || z >= def_Nz-1u) return;

  // Get indices of 8 cube corners
  const uxx x0 = (uxx)x, xp = (uxx)(x+1u);
  const uxx y0 = (uxx)(y*def_Nx), yp = (uxx)((y+1u)*def_Nx);
  const uxx z0 = (uxx)(z)*(uxx)(def_Ny*def_Nx), zp = (uxx)(z+1u)*(uxx)(def_Ny*def_Nx);

#ifdef SURFACE_EXPORT_IGNORE_BOUNDARIES
  // Check if any of the 8 cube corners is a solid boundary cell
  const uxx corners[8] = {
    x0+y0+z0, xp+y0+z0, xp+y0+zp, x0+y0+zp,
    x0+yp+z0, xp+yp+z0, xp+yp+zp, x0+yp+zp
  };
  for(uint i = 0u; i < 8u; i++) {
    if(flags[corners[i]] & TYPE_S) return; // Skip this cube if any corner touches a solid boundary
  }
#endif

  float v[8];
  v[0] = phi[x0+y0+z0]; v[1] = phi[xp+y0+z0];
  v[2] = phi[xp+y0+zp]; v[3] = phi[x0+y0+zp];
  v[4] = phi[x0+yp+z0]; v[5] = phi[xp+yp+z0];
  v[6] = phi[xp+yp+zp]; v[7] = phi[x0+yp+zp];

  // Get triangles from marching cubes
  float3 triangles[15];
  const uint tn = marching_cubes(v, 0.5f, triangles);

  if(tn > 0u) {
    // Atomically reserve space in output buffer
    const uint base_index = atomic_add(triangle_count, tn);

    // Check buffer overflow
    if(base_index + tn > max_triangles) return;

    // Write triangles to global memory
    const float3 offset = (float3)((float)x, (float)y, (float)z);
    for(uint i = 0u; i < tn; i++) {
      const uint vertex_index = (base_index + i) * 9u;
      const float3 p0 = triangles[3u*i] + offset;
      const float3 p1 = triangles[3u*i+1u] + offset;
      const float3 p2 = triangles[3u*i+2u] + offset;

      // Write triangle vertices (9 floats per triangle)
      vertices[vertex_index+0u] = p0.x;
      vertices[vertex_index+1u] = p0.y;
      vertices[vertex_index+2u] = p0.z;
      vertices[vertex_index+3u] = p1.x;
      vertices[vertex_index+4u] = p1.y;
      vertices[vertex_index+5u] = p1.z;
      vertices[vertex_index+6u] = p2.x;
      vertices[vertex_index+7u] = p2.y;
      vertices[vertex_index+8u] = p2.z;
    }
  }
}
#endif // SURFACE_EXPORT

[nhemsley]

OK, I have extracted this functionality into https://github.com/nhemsley/FluidX3D/tree/surface-export

This has one commit:

6332ed5

This was admittedly written with alot of help by claude, and then re-extracted. This is therefore a proof of concept. If you have any changes you wish to make, I can help out to refine it.

Notes:

setup.cpp. I have left the beach example uncommented.

Marching Cubes:

Marching cubes code is only defined, when GRAPHICS is defined. This probably should be gated in i.e. MARCHING_CUBES (or maybe SURFACE ?). I have left this as is at the moment.

surface_export.cpp:

SurfaceExportConfig, This may not belong in FluidX3D proper, I can remove this if you would prefer.
These other stl functions may want to go into utilities.cpp.

Probably some other but I will go over the commit again at a later point, and mention anything that comes to mind.

[nhemsley]

Here is my defines, to enable the settings properly, this will export to FluidX3D/export/beach_waves. notably:

#define SURFACE_EXPORT // enables GPU kernels for extracting surface mesh via marching cubes and exporting to STL files; requires SURFACE & GRAPHICS

#pragma once



//#define D2Q9 // choose D2Q9 velocity set for 2D; allocates 53 (FP32) or 35 (FP16) Bytes/cell
//#define D3Q15 // choose D3Q15 velocity set for 3D; allocates 77 (FP32) or 47 (FP16) Bytes/cell
#define D3Q19 // choose D3Q19 velocity set for 3D; allocates 93 (FP32) or 55 (FP16) Bytes/cell; (default)
//#define D3Q27 // choose D3Q27 velocity set for 3D; allocates 125 (FP32) or 71 (FP16) Bytes/cell

#define SRT // choose single-relaxation-time LBM collision operator; (default)
//#define TRT // choose two-relaxation-time LBM collision operator

#define FP16S // optional for 2x speedup and 2x VRAM footprint reduction: compress LBM DDFs to range-shifted IEEE-754 FP16; number conversion is done in hardware; all arithmetic is still done in FP32
//#define FP16C // optional for 2x speedup and 2x VRAM footprint reduction: compress LBM DDFs to more accurate custom FP16C format; number conversion is emulated in software; all arithmetic is still done in FP32

// #define BENCHMARK // disable all extensions and setups and run benchmark setup instead

#define VOLUME_FORCE // enables global force per volume in one direction (equivalent to a pressure gradient); specified in the LBM class constructor; the force can be changed on-the-fly between time steps at no performance cost
//#define FORCE_FIELD // enables computing the forces on solid boundaries with lbm.update_force_field(); and enables setting the force for each lattice point independently (enable VOLUME_FORCE too); allocates an extra 12 Bytes/cell
#define EQUILIBRIUM_BOUNDARIES // enables fixing the velocity/density by marking cells with TYPE_E; can be used for inflow/outflow; does not reflect shock waves
//#define MOVING_BOUNDARIES // enables moving solids: set solid cells to TYPE_S and set their velocity u unequal to zero
#define SURFACE // enables free surface LBM: mark fluid cells with TYPE_F; at initialization the TYPE_I interface and TYPE_G gas domains will automatically be completed; allocates an extra 12 Bytes/cell
#define SURFACE_EXPORT // enables GPU kernels for extracting surface mesh via marching cubes and exporting to STL files; requires SURFACE & GRAPHICS
#define SURFACE_EXPORT_IGNORE_BOUNDARIES // optional: skip exporting triangles where surface touches solid boundaries (TYPE_S); reduces triangle count significantly for enclosed domains
//#define TEMPERATURE // enables temperature extension; set fixed-temperature cells with TYPE_T (similar to EQUILIBRIUM_BOUNDARIES); allocates an extra 32 (FP32) or 18 (FP16) Bytes/cell
//#define SUBGRID // enables Smagorinsky-Lilly subgrid turbulence LES model to keep simulations with very large Reynolds number stable
//#define PARTICLES // enables particles with immersed-boundary method (for 2-way coupling also activate VOLUME_FORCE and FORCE_FIELD; only supported in single-GPU)

//#define INTERACTIVE_GRAPHICS // enable interactive graphics; start/pause the simulation by pressing P; either Windows or Linux X11 desktop must be available; on Linux: change to "compile on Linux with X11" command in make.sh
//#define INTERACTIVE_GRAPHICS_ASCII // enable interactive graphics in ASCII mode the console; start/pause the simulation by pressing P
#define GRAPHICS // run FluidX3D in the console, but still enable graphics functionality for writing rendered frames to the hard drive

#define GRAPHICS_FRAME_WIDTH 1920 // set frame width if only GRAPHICS is enabled
#define GRAPHICS_FRAME_HEIGHT 1080 // set frame height if only GRAPHICS is enabled
#define GRAPHICS_BACKGROUND_COLOR 0x000000 // set background color; black background (default) = 0x000000, white background = 0xFFFFFF
#define GRAPHICS_U_MAX 0.18f // maximum velocity for velocity coloring in units of LBM lattice speed of sound (c=1/sqrt(3)) (default: 0.18f)
#define GRAPHICS_RHO_DELTA 0.001f // coloring range for density rho will be [1.0f-GRAPHICS_RHO_DELTA, 1.0f+GRAPHICS_RHO_DELTA] (default: 0.001f)
#define GRAPHICS_T_DELTA 1.0f // coloring range for temperature T will be [1.0f-GRAPHICS_T_DELTA, 1.0f+GRAPHICS_T_DELTA] (default: 1.0f)
#define GRAPHICS_F_MAX 0.001f // maximum force in LBM units for visualization of forces on solid boundaries if VOLUME_FORCE is enabled and lbm.update_force_field(); is called (default: 0.001f)
#define GRAPHICS_Q_CRITERION 0.0001f // Q-criterion value for Q-criterion isosurface visualization (default: 0.0001f)
#define GRAPHICS_STREAMLINE_SPARSE 8 // set how many streamlines there are every x lattice points
#define GRAPHICS_STREAMLINE_LENGTH 128 // set maximum length of streamlines
#define GRAPHICS_RAYTRACING_TRANSMITTANCE 0.25f // transmitted light fraction in raytracing graphics ("0.25f" = 1/4 of light is transmitted and 3/4 is absorbed along longest box side length, "1.0f" = no absorption)
#define GRAPHICS_RAYTRACING_COLOR 0x005F7F // absorption color of fluid in raytracing graphics

//#define GRAPHICS_TRANSPARENCY 0.7f // optional: comment/uncomment this line to disable/enable semi-transparent rendering (looks better but reduces framerate), number represents transparency (equal to 1-opacity) (default: 0.7f)



// #############################################################################################################

#define TYPE_S 0b00000001 // (stationary or moving) solid boundary
#define TYPE_E 0b00000010 // equilibrium boundary (inflow/outflow)
#define TYPE_T 0b00000100 // temperature boundary
#define TYPE_F 0b00001000 // fluid
#define TYPE_I 0b00010000 // interface
#define TYPE_G 0b00100000 // gas
#define TYPE_X 0b01000000 // reserved type X
#define TYPE_Y 0b10000000 // reserved type Y

#define VIS_FLAG_LATTICE  0b00000001 // lbm.graphics.visualization_modes = VIS_...|VIS_...|VIS_...;
#define VIS_FLAG_SURFACE  0b00000010
#define VIS_FIELD         0b00000100
#define VIS_STREAMLINES   0b00001000
#define VIS_Q_CRITERION   0b00010000
#define VIS_PHI_RASTERIZE 0b00100000
#define VIS_PHI_RAYTRACE  0b01000000
#define VIS_PARTICLES     0b10000000

#if defined(FP16S) || defined(FP16C)
#define fpxx ushort
#else // FP32
#define fpxx float
#endif // FP32

#ifdef BENCHMARK
#undef UPDATE_FIELDS
#undef VOLUME_FORCE
#undef FORCE_FIELD
#undef MOVING_BOUNDARIES
#undef EQUILIBRIUM_BOUNDARIES
#undef SURFACE
#undef TEMPERATURE
#undef SUBGRID
#undef PARTICLES
#undef INTERACTIVE_GRAPHICS
#undef INTERACTIVE_GRAPHICS_ASCII
#undef GRAPHICS
#endif // BENCHMARK

#ifdef SURFACE // (rho, u) need to be updated exactly every LBM step
#define UPDATE_FIELDS // update (rho, u, T) in every LBM step
#endif // SURFACE

#ifdef TEMPERATURE
#define VOLUME_FORCE
#endif // TEMPERATURE

#ifdef PARTICLES // (rho, u) need to be updated exactly every LBM step
#define UPDATE_FIELDS // update (rho, u, T) in every LBM step
#endif // PARTICLES

#if defined(INTERACTIVE_GRAPHICS) || defined(INTERACTIVE_GRAPHICS_ASCII)
#define GRAPHICS
#define UPDATE_FIELDS // to prevent flickering artifacts in interactive graphics
#endif // INTERACTIVE_GRAPHICS || INTERACTIVE_GRAPHICS_ASCII