Equilibrium Boundaries with Fluids

[C-Dub2022]

I can't for the life of me seem to get equilibrium boundaries to work with surfaces. I just seem to end up with the surface level at the front of the volume slowly dropping to nothing. Any ideas? I'm probably doing something very silly. This is my initial conditions. I'm probably doing something very silly.....

if (z < H && lbm.flags[n] != TYPE_S) {
	lbm.flags[n] = TYPE_F;
	lbm.rho[n] = units.rho_hydrostatic(units.f(si_rho, si_g), (float)z, (float)H);
	lbm.u.y[n] = -lbm_u;
}

if (x == 0u || x == Nx - 1u || z == Nz - 1u || y == 0u || y == Ny - 1u) lbm.flags[n] = TYPE_E; // all non periodic

if (y == Ny - 1u) {
	lbm.rho[n] = 1.0f * units.rho_hydrostatic(units.f(si_rho, si_g), (float)z, (float)H);
	lbm.u.y[n] = -lbm_u;
}

if (z == 0u) lbm.flags[n] = TYPE_S;

[ProjectPhysX]

Hi @C-Dub2022, the one small mistake is setting lbm.flags[n] = TYPE_E;. This clears the previously set TYPE_F flag in the boundary cells. Instead, set lbm.flags[n] |= TYPE_E;, note the |; this sets the TYPE_E flag additionally to the previously set TYPE_F flag. Please let me know if this fixes the problem in your setup.

Here is a sample setup for testing:

void main_setup() { // hydrodynamics of a cow; required extensions in defines.hpp: FP16S, EQUILIBRIUM_BOUNDARIES, SURFACE, SUBGRID, INTERACTIVE_GRAPHICS or GRAPHICS
	// ################################################################## define simulation box size, viscosity and volume force ###################################################################
	const uint3 lbm_N = resolution(float3(1.0f, 2.0f, 1.0f), 1000u); // input: simulation box aspect ratio and VRAM occupation in MB, output: grid resolution
	const float si_u = 1.0f;
	const float si_length = 2.4f;
	const float si_T = 10.0f;
	const float si_nu=1.0E-6f, si_rho=1000.0f;
	const float si_g = 9.81f;
	const float lbm_length = 0.65f*(float)lbm_N.y;
	const float lbm_h = 0.5f*(float)lbm_N.z;
	const float lbm_u = 0.1f;
	units.set_m_kg_s(lbm_length, lbm_u, 1.0f, si_length, si_u, si_rho);
	print_info("Re = "+to_string(to_uint(units.si_Re(si_length, si_u, si_nu))));
	const float lbm_f = units.f(si_rho, si_g);
	LBM lbm(lbm_N, units.nu(si_nu), 0.0f, 0.0f, -lbm_f);
	// ###################################################################################### define geometry ######################################################################################
	const float3x3 rotation = float3x3(float3(1, 0, 0), radians(180.0f))*float3x3(float3(0, 0, 1), radians(180.0f));
	Mesh* mesh = read_stl(get_exe_path()+"../stl/Cow_t.stl", lbm.size(), lbm.center(), rotation, lbm_length); // https://www.thingiverse.com/thing:182114/files
	mesh->translate(float3(0.0f, 1.0f-mesh->pmin.y+0.1f*lbm_length, 1.0f-mesh->pmin.z)); // move mesh forward a bit and to simulation box bottom, keep in mind 1 cell thick box boundaries
	lbm.voxelize_mesh_on_device(mesh);
	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);
		if(z==0u) lbm.flags[n] = TYPE_S; // solid floor
		if(lbm.flags[n]!=TYPE_S) {
			lbm.u.y[n] = lbm_u; // initialize y-velocity everywhere except in solid cells
			if(z<=lbm_h) {
				if(lbm.flags[n]!=TYPE_S) lbm.flags[n] = TYPE_F;
				lbm.rho[n] = units.rho_hydrostatic(lbm_f, (float)z, lbm_h);
			}
		}
		if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==Nz-1u) lbm.flags[n] |= TYPE_E; // all other simulation box boundaries are inflow/outflow
	}); // ####################################################################### run simulation, export images and data ##########################################################################
	lbm.graphics.visualization_modes = VIS_FLAG_LATTICE|VIS_FLAG_SURFACE|VIS_Q_CRITERION|VIS_PHI_RASTERIZE;
#if defined(GRAPHICS) && !defined(INTERACTIVE_GRAPHICS)
	lbm.graphics.set_camera_centered(-40.0f, 20.0f, 78.0f, 1.25f);
	lbm.run(0u); // initialize simulation
	while(lbm.get_t()<=units.t(si_T)) { // main simulation loop
		if(lbm.graphics.next_frame(units.t(si_T), 10.0f)) lbm.graphics.write_frame();
		lbm.run(1u);
	}
#else // GRAPHICS && !INTERACTIVE_GRAPHICS
	lbm.run();
#endif // GRAPHICS && !INTERACTIVE_GRAPHICS
} /**/