ProjectPhysX
Question about Solving for stability
Hello,
I'm running a simulation with an operating Reynolds number set to approximately 2×10^7. However, I quickly noticed the rapid growth of white, crystal-like formations on the bottom wall and around the structure. Based on your earlier description, it appears that these artifacts are due to simulation instability.
For reference, my structure has a characteristic length of 30 m, a velocity of 10 m/s, an air viscosity of 1.48×10⁻⁵, and I've set lbm_u to 0.05f. Could you please advise if you have any suggestions to resolve this issue? I've been struggling with it for quite some time.
Thank you very much for your help.
Dear Moritz,
Thank you for your quick response.
I enabled the SUBGRID extension. I attached my project. Could u give some advices for me? Thank you very much.
Best regards, Pengfei
发件人: Dr. Moritz Lehmann @.> 发送时间: 2025年2月8日 22:44 收件人: ProjectPhysX/FluidX3D @.> 抄送: LIN Pengfei @.>; Author @.> 主题: Re: [ProjectPhysX/FluidX3D] Question about Solving for stability (Issue #268)
Hi @ProjectPhysXhttps://github.com/ProjectPhysX,
have you enabled the SUBGRIDhttps://github.com/ProjectPhysX/FluidX3D/blob/master/src/defines.hpp#L24 extension?
Kind regards, Moritz
― Reply to this email directly, view it on GitHubhttps://github.com/ProjectPhysX/FluidX3D/issues/268#issuecomment-2645748263, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AZA6ARB3BD3UIKYZGVOKJ532OYJ5BAVCNFSM6AAAAABWXRH7U2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNBVG42DQMRWGM. You are receiving this because you authored the thread.Message ID: @.***>
#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.calculate_force_on_boundaries(); 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 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.20f // 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.01f // coloring range for density rho will be [1.0f-GRAPHICS_RHO_DELTA, 1.0f+GRAPHICS_RHO_DELTA] (default: 0.01f) #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.calculate_force_on_boundaries(); 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
void main_setup() { // required extensions in defines.hpp: FP16C, FORCE_FIELD, EQUILIBRIUM_BOUNDARIES, SUBGRID, optionally INTERACTIVE_GRAPHICS
const uint memory = 5500u;
const float lbm_u = 0.05f;
const float si_u = 10.0f;
const float si_R_height = 0.1f;
const float si_nu = 1.48E-5f, si_rho = 1.225f;
const float si_width = 30.0f, si_height = 180.0f, si_length = 30.0f;
const float si_A = si_width * si_height;
const float si_t = 60.0f;
int n_jpg = 0;
const uint3 lbm_N = resolution(float3(1.0f, 2.4f, 1.0f), memory);
units.set_m_kg_s((float)lbm_N.y, lbm_u, 1.225f, 7.37f, si_u, si_rho);
const float lbm_nu = units.nu(si_nu);
const ulong lbm_T = units.t(si_t);
print_info("lbm_nu=" + to_string((lbm_nu)));
print_info("lbm_T=" + to_string((lbm_T)));
print_info("lbm_N.x=" + to_string(((float)lbm_N.x)));
print_info("lbm_N.y=" + to_string(((float)lbm_N.y)));
print_info("lbm_N.z=" + to_string(((float)lbm_N.z)));
print_info("Re = " + to_string(to_uint(units.si_Re(si_width, si_u, si_nu))));
LBM lbm(lbm_N, lbm_nu);
const float lbm_width = units.x(si_width);
const float lbm_height = units.x(si_height);
const float lbm_R_height = units.x(si_R_height);
Mesh* mesh = read_stl(get_exe_path() + "stl/30000.stl", lbm.size(), lbm.center(), float3x3(float3(0, 0, 1), radians(0.0f)), lbm_height);
mesh->translate(float3(0.0f, (5.50f * lbm_width - 0.5f * (float)lbm_N.y), 0.5f * lbm_height - lbm.center().z));
lbm.voxelize_mesh_on_device(mesh, TYPE_S | TYPE_X);
const string path = get_exe_path() + "../result/add-tur-01/";
write_file(path + "Cd.dat", "# t\tCd\n");
parallel_for(lbm.get_N(), [&](ulong n) {
uint x = 0u, y = 0u, z = 0u; lbm.coordinates(n, x, y, z);
const float wind_speed = lbm_u * std::pow(z / lbm_R_height, 0.22);
if (lbm.flags[n] != TYPE_S && lbm.flags[n] != (TYPE_S | TYPE_X)) {
lbm.u.y[n] = wind_speed;
lbm.u.x[n] = 0.0f;
lbm.u.z[n] = 0.0f;
};
if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == Nz - 1u) lbm.flags[n] = TYPE_E;
if (z == 0u) lbm.flags[n] = TYPE_S;});
lbm.graphics.visualization_modes = VIS_FLAG_SURFACE | VIS_Q_CRITERION;
lbm.graphics.set_camera_centered(300.0f, 400.0f, 100.0f, 20.0f);
lbm.run(0u, lbm_T);
lbm.write_status(path);
while (lbm.get_t() <= lbm_T) {
if (lbm.graphics.next_frame(lbm_T, 100.0f)) {
lbm.calculate_force_on_boundaries();
lbm.F.read_from_device();
//lbm.u.write_device_to_vtk(); // velocity
//lbm.F.write_device_to_vtk(); // force, only for FORCE_FIELD extension
lbm.voxelize_mesh_on_device(mesh, TYPE_S | TYPE_X);
const float3 lbm_force = lbm.calculate_force_on_object(TYPE_S | TYPE_X);
float si_u_B_height = units.si_u(lbm_u * std::pow(lbm_height / lbm_R_height, 0.22));
const float Cd = units.si_F(lbm_force.y) / (0.5f * si_rho * sq(si_u_B_height) * si_A);
const float Cl = units.si_F(lbm_force.x) / (0.5f * si_rho * sq(si_u_B_height) * si_A);
print_info("\r" + to_string(Cd, 10u) + "Cd" + to_string(clock.stop(), 3u) + "Time");
print_info("\r" + to_string(Cl, 10u) + "Cl" + to_string(clock.stop(), 3u) + "Time");
write_line(path + "Cd.dat", to_string(lbm.get_t()) + "\t" + to_string(Cd, 10u) + "\n");
lbm.graphics.write_frame_png(get_exe_path() + "jpg/square-30-add-tur-01/");
n_jpg++;
}
lbm.run(1u, lbm_T);
}
lbm.run();
}
Hi Moritz,
Meanwhile, in my subsequent data processing, I discovered that when the inlet velocity impacts the surface of the structure, the wind speed bounces back in the opposite direction, and its magnitude is relatively high (for example, when the inlet is 20 m/s, the rebound speed exceeds 10 m/s). This discrepancy from actual behavior is too significant! Could you offer some suggestions? I would also be very grateful for any advice you might have regarding my previous question.
Best regards, Pengfei