D3Q7 thermal stability at τ_T → 0.5 - any plans for TRT on the temperature lattice?

Run the setup.cpp reproducer with the listed TEMPERATURE, FP16S, and EQUILIBRIUM_BOUNDARIES defines, and record the D3Q7 overshoot and NaN behavior near τ_T≈0.5. Read defines.hpp and the temperature-related implementation, then define the accepted stabilization scope and verify that the case no longer diverges without relying on the per-step clamp.

Dear Dr. Lehmann,
First, thank you for FluidX3D — the performance is remarkable, and we've successfully used it for pressure-drop and flow-distribution analysis of microchannel water cold plates.
I'd like to ask about the temperature extension. Our case: water (Pr ≈ 3 at film temperature) in 0.1 mm channels, D3Q19 SRT + FP16S, dx = 6.25 µm, with TYPE_T Dirichlet cells on the wetted surface. Realistic velocity/pressure scaling forces lbm_u ≈ 0.008, which puts the thermal relaxation at τ_T ≈ 0.504 (w_T ≈ 1.985). The hydrodynamic solution converges cleanly (mass-balanced, dP plausible), but the D3Q7 scalar field rings strongly at the sharp Dirichlet boundaries (±100 K over/undershoot) and eventually diverges to NaN. A per-step clamp on T stabilizes it, but that's obviously a workaround. One side observation: once T is NaN, the Boussinesq force term propagates NaN into the flow even with beta = 0 (0 × NaN), so the hydrodynamic field is lost as well.
Is a TRT/two-relaxation option (or another stabilization) for the D3Q7 temperature lattice something you've considered or would accept a contribution toward? Or would you recommend a different parameterization for advection-dominated thermal cases like this one? I have a compact reproducible setup if useful.
Best regards,
Vu

setup.cpp

void main_setup() { // v3 (sinusoidal fin array, periodic-strip, split-flow along y); required extensions in defines.hpp: FP16S, EQUILIBRIUM_BOUNDARIES, TEMPERATURE, INTERACTIVE_GRAPHICS
	// ============================================================= geometry (measured from meshed_fin_array.stl) =============================================================
	// plate 50.5 x 50.3 x 2.0 mm; fin field 28.0 x 27.0 mm centered; fins 5.0 mm tall; continuous sinusoidal walls along y, periodic pitch 0.301 mm along x.
	// Domain = strip of 1 fin pitch in x with periodic x boundaries (exact for the array interior; voxelization verified watertight and
	// fully connected inlet->outlets). y is cropped to the fin field + small plenums; outlets sit at the channel ends as designed.
	// v3 changes vs v2: (1) lbm_u lowered 0.10 -> 0.025: the ~1 kPa channel pressure drop in lattice units scales with 1/lbm_u^2;
	// at lbm_u=0.10 it pushed rho to ~1.2-1.7 -> FP16S DDF overflow -> NaN ("N/A") once the channels filled and pressure reflected.
	// (2) finer grid (32 cells/pitch) to keep tau_T off the 0.5 limit despite the slower lattice velocity. (3) inlet sits directly
	// on the fin tips (design intent: no cool water wasted above the array). (4) codex-style multi-line console diagnostics.
	const float si_pitch        = 0.2E-3f; // [m] fin pitch along x, measured from straight_fins.stl (0.1 fin + 0.1 channel); must equal the loaded STL's repeat distance
	const uint  cells_per_pitch = 32u;       // -> dx = 6.25 um, ~16 cells across each 0.1 mm channel (also raises tau_T at low lbm_u)
	const uint  pitches         = 1u;        // strip width in pitches (x periodic); increase when more VRAM is available
	const float si_fin_len_x    = 50.8E-3f;  // [m] fin-field extent along x (straight fins span the full plate; sets the strip's share of the total flow)
	const float si_fin_len_y    = 27.0E-3f;  // [m] fin length along y = channel length
	const float si_fin_height   = 2.45E-3f;  // [m] measured from straight_fins.stl
	const float si_base_height  = 2.7E-3f;   // [m] measured from straight_fins.stl; plate bottom on the adiabatic domain floor
	const float si_outlet_len   = 0.5E-3f;   // [m] thin outlet plenum right at each channel end (y faces = outlets)
	const float si_dx           = si_pitch/(float)cells_per_pitch; // [m] lattice spacing
	const uint Nx = pitches*cells_per_pitch;
	const uint Ny = ((uint)((si_fin_len_y+2.0f*si_outlet_len)/si_dx)+1u)&~1u;
	const uint Nz = (((uint)ceil((si_base_height+si_fin_height)/si_dx)+3u)&~1u); // fins reach to ~2 cells below the lid; inlet row sits directly on the fin tips

	// ============================================================= SI parameters =============================================================
	const float si_T_in   = 25.0f+273.15f;   // [K] coolant inlet temperature
	const float si_T_wall = 90.0f+273.15f;   // [K] isothermal wetted-surface proxy. WARNING: no solid conduction in FluidX3D -> no fin efficiency -> Q is an UPPER BOUND.
	const float si_T_film = 0.5f*(si_T_in+si_T_wall); // ~57.5 degC
	// water properties at the film temperature:
	const float si_rho   = 984.3f;           // [kg/m^3]
	const float si_cp    = 4184.0f;          // [J/(kg*K)]
	const float si_nu    = 0.488E-6f;        // [m^2/s]
	const float si_Pr    = 3.08f;
	const float si_alpha = si_nu/si_Pr;      // [m^2/s]
	const float si_beta  = 0.0f;             // forced convection; buoyancy off
	const float si_Q_total = 2.55f/60000.0f; // [m^3/s] = 2.55 L/min through the whole heat sink
	const float si_strip_w = (float)pitches*si_pitch;            // [m]
	const float si_Q       = si_Q_total*si_strip_w/si_fin_len_x; // [m^3/s] strip share
	const float si_inlet_w = 8.0E-3f;        // [m] center manifold slot width along y (slot spans all of x = covers every channel)
	const float si_u_in    = si_Q/(si_inlet_w*si_strip_w);       // [m/s] downward inlet velocity (~0.19 m/s)

	// ============================================================= LBM units =============================================================
	// pressure unit = si_rho*(si_u_in/lbm_u)^2: lbm_u=0.025 -> ~57 kPa per LBM pressure unit -> the ~1-2 kPa physical dP stays at
	// rho deviation <~0.06 (weakly compressible, safe for FP16S). Trade-off: tau_T = 0.5+2*alpha ~ 0.504 (watch the overshoot warning).
	// lbm_u sets the pressure headroom: p_unit = si_rho*(si_u_in/lbm_u)^2 must be >> the real pump head this geometry needs.
	// Straight 2.45 mm fins need ~4 kPa to pass 2.55 L/min (U_ch ~ 0.68 m/s) -> lbm_u = 0.008 gives p_unit ~ 170 kPa (drho ~ 0.07).
	const float lbm_u = 0.008f;
	units.set_m_kg_s_K((float)Nx, lbm_u, 1.0f, 1.0f, (float)Nx*si_dx, si_u_in, si_rho, si_T_film);
	const float lbm_nu     = units.nu(si_nu);
	const float lbm_alpha  = units.alpha(si_alpha);
	const float lbm_beta   = units.beta(si_beta);
	const float lbm_T_in   = units.T(si_T_in);
	const float lbm_T_wall = units.T(si_T_wall);

	LBM lbm(Nx, Ny, Nz, lbm_nu, 0.0f, 0.0f, 0.0f, 0.0f, lbm_alpha, lbm_beta);

	// ============================================================= voxelize geometry =============================================================
	Mesh* heatsink = read_stl(get_exe_path()+"../stl/straight_fins.stl");
	heatsink->scale(0.001f);        // mm -> m
	heatsink->scale(units.x(1.0f)); // m -> lattice cells
	const float3 mesh_center = 0.5f*(heatsink->pmin+heatsink->pmax);
	heatsink->translate(float3(lbm.center().x-mesh_center.x, lbm.center().y-mesh_center.y, 0.0f-heatsink->pmin.z)); // plate bottom on the domain floor
	lbm.voxelize_mesh_on_device(heatsink, TYPE_S); // z-ray voxelization (forced in lbm.cpp); verified watertight for this mesh
	lbm.flags.read_from_device();

	const uint hot_top_z = (uint)ceil((si_base_height+0.90f*si_fin_height)/si_dx); // top 10% of fin height stays adiabatic: real fin tips are the coldest part of the fin (finite fin efficiency), and removing the 90 degC Dirichlet cells from the inlet jet-shear zone eliminates the main D3Q7 overshoot source
	const float yc = lbm.center().y;
	const float lbm_inlet_half_w = 0.5f*si_inlet_w/si_dx;
	const float ramp_w = 0.5E-3f/si_dx; // smooth the inlet velocity over 0.5 mm at the slot edges (removes the sharp shear corner at the slot edge that drives local overshoot)
	double ramp_sum = 0.0; ulong ramp_n = 0ull;
	for(uint y=0u; y<Ny; y++) {
		const float sr = (lbm_inlet_half_w-fabs((float)y+0.5f-yc))/ramp_w;
		if(sr>0.0f) { const float r = fmin(sr, 1.0f); ramp_sum += (double)(r*r*(3.0f-2.0f*r)); ramp_n++; }
	}
	const float ramp_norm = ramp_n ? (float)((double)ramp_n/ramp_sum) : 1.0f; // rescale so the total flow rate stays exact

	// hot wall: TYPE_T on FLUID cells adjacent to the solid surface = FluidX3D's thermal Dirichlet BC (solids themselves are adiabatic bounce-back)
	vector<uchar> hot_wall(lbm.get_N(), 0u);
	parallel_for(lbm.get_N(), [&](ulong n) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
		if((lbm.flags[n]&TYPE_S)==TYPE_S || y==0u || y==Ny-1u || z==0u || z==Nz-1u) return;
		if(z>hot_top_z) return;
		const uint xm = (x+Nx-1u)%Nx, xp = (x+1u)%Nx; // x is periodic
		const bool adj =
			((lbm.flags[lbm.index(xm, y, z)]&TYPE_S)==TYPE_S) || ((lbm.flags[lbm.index(xp, y, z)]&TYPE_S)==TYPE_S) ||
			((lbm.flags[lbm.index(x, y-1u, z)]&TYPE_S)==TYPE_S) || ((lbm.flags[lbm.index(x, y+1u, z)]&TYPE_S)==TYPE_S) ||
			((lbm.flags[lbm.index(x, y, z-1u)]&TYPE_S)==TYPE_S) || ((lbm.flags[lbm.index(x, y, z+1u)]&TYPE_S)==TYPE_S);
		hot_wall[n] = adj ? 1u : 0u;
	});

	parallel_for(lbm.get_N(), [&](ulong n) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
		lbm.T[n] = lbm_T_in; // cold start
		if((lbm.flags[n]&TYPE_S)==TYPE_S) return;
		const bool in_slot = fabs((float)y+0.5f-yc)<=lbm_inlet_half_w;
		if(z==Nz-1u) {
			if(in_slot) { // manifold slot directly on the fin tips (design intent); smoothstep taper at the slot edges, normalized to keep Q exact
				const float sr = fmin((lbm_inlet_half_w-fabs((float)y+0.5f-yc))/ramp_w, 1.0f);
				lbm.flags[n] = TYPE_E|TYPE_T;
				lbm.u.z[n] = -lbm_u*ramp_norm*sr*sr*(3.0f-2.0f*sr);
			}
			else lbm.flags[n] = TYPE_S; // closed lid elsewhere
		} else if(z==0u) {
			lbm.flags[n] = TYPE_S; // adiabatic floor
		} else if(y==0u || y==Ny-1u) {
			lbm.flags[n] = TYPE_E; // pressure outlets (rho=1, u=0) at the channel ends; temperature floats (soft outflow)
		} else if(hot_wall[n]) {
			lbm.flags[n] = TYPE_T; lbm.T[n] = lbm_T_wall; // isothermal wetted surface
		}
		// x faces: no flags -> periodic
	});

	ulong solid_count=0ull, hot_count=0ull, inlet_count=0ull, outlet_count=0ull;
	for(ulong n=0ull; n<lbm.get_N(); n++) {
		uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
		solid_count  += (lbm.flags[n]&TYPE_S)==TYPE_S;
		hot_count    += (lbm.flags[n]&(TYPE_S|TYPE_E|TYPE_T))==TYPE_T;
		inlet_count  += ((lbm.flags[n]&(TYPE_E|TYPE_T))==(TYPE_E|TYPE_T))&&z==Nz-1u;
		outlet_count += ((lbm.flags[n]&(TYPE_S|TYPE_E|TYPE_T))==TYPE_E)&&(y==0u||y==Ny-1u);
	}

	// ============================================================= info =============================================================
	const float si_u_ch_est = si_Q_total/(si_fin_len_x*si_fin_height*(1.0f-0.1E-3f/si_pitch));
	print_info("=========================== fin-array strip v3 ===========================");
	print_info("Domain        : "+to_string(Nx)+" x "+to_string(Ny)+" x "+to_string(Nz)+" cells = "+to_string((float)(lbm.get_N()/1000000ull), 0u)+"M ("+to_string(pitches)+" pitch, x periodic), dx = "+to_string(1E6f*si_dx, 2)+" um, ~"+to_string(0.1E-3f/si_dx, 1)+" cells across a 0.1 mm gap");
	print_info("Cells         : solid = "+to_string(solid_count)+", hot-wall layer = "+to_string(hot_count)+", inlet = "+to_string(inlet_count)+", outlets = "+to_string(outlet_count));
	print_info("Flow          : Q_total = "+to_string(si_Q_total*60000.0f, 2)+" L/min, strip share = "+to_string(si_Q*1E9f, 1)+" uL/s, U_in = "+to_string(si_u_in, 3)+" m/s, U_ch ~ "+to_string(si_u_ch_est, 3)+" m/s, Re_ch ~ "+to_string(to_uint(si_u_ch_est*2.0f*0.1E-3f/si_nu))+" (laminar)");
	print_info("Water @ 57.5C : nu = 0.488E-6 m^2/s, Pr = "+to_string(si_Pr, 2)+" | T_in = "+to_string(si_T_in-273.15f, 1)+" degC, T_wall = "+to_string(si_T_wall-273.15f, 1)+" degC (isothermal proxy -> Q = upper bound)");
	print_info("Numerics      : lbm_u = "+to_string(lbm_u, 3)+", tau = "+to_string(3.0f*lbm_nu+0.5f, 5)+", tau_T = "+to_string(2.0f*lbm_alpha+0.5f, 5)+" (w_T = "+to_string(1.0f/(2.0f*lbm_alpha+0.5f), 4)+" / limit 2.0), p_unit = "+to_string(units.si_p(1.0f)/1000.0f, 1)+" kPa"
#ifdef TEMPERATURE_LIMIT
		+" | T-limiter ON ["+to_string(units.si_T(1.0f-TEMPERATURE_LIMIT)-273.15f, 1)+", "+to_string(units.si_T(1.0f+TEMPERATURE_LIMIT)-273.15f, 1)+"] degC"
#else
		+" | T-limiter OFF"
#endif
	);
	print_info("==========================================================================");

	// ============================================================= run =============================================================
	lbm.graphics.visualization_modes = VIS_FIELD|VIS_STREAMLINES;
	lbm.graphics.field_mode = 2; // temperature; press Z to cycle velocity/pressure/temperature
	lbm.graphics.slice_mode = 1;
	lbm.graphics.set_camera_centered(0.0f, 20.0f, 60.0f, 1.2f);

	const ulong dt = 4000ull;       // diagnostics every dt steps
	const ulong t_min = 200000ull;  // no convergence verdict before this (several channel transit times)
	const ulong t_max = 3000000ull; // hard stop
	const bool export_vtk = false;  // set true to export rho/u/T/flags every vtk_every steps
	const ulong vtk_every = 200000ull;
	ulong next_vtk = vtk_every;
	const string vtk_path = get_exe_path()+"export/vtk/";

	double Tout2=0.0, Tout1=0.0, drho_in=0.0; uint conv_samples=0u;
	while(lbm.get_t()<t_max) {
		lbm.run(dt);
		lbm.rho.read_from_device(); lbm.u.read_from_device(); lbm.T.read_from_device();

		double sum_out_uT=0.0, sum_out_u=0.0, sum_in_u=0.0, sum_in_p=0.0, sum_out_p=0.0, sum_near_T=0.0, sum_in_uz_real=0.0, max_u=0.0, max_drho=0.0;
		ulong p_in_n=0ull, p_out_n=0ull, near_n=0ull;
		float min_Tf=max_float, max_Tf=-max_float, min_Tb=max_float, max_Tb=-max_float;
		ulong nTf=0ull, nTb=0ull;
		uint mux=0u, muy=0u, muz=0u, mTx=0u, mTy=0u, mTz=0u; uchar muf=0u, mTf2=0u; float max_T_all=-max_float;
		for(ulong n=0ull; n<lbm.get_N(); n++) {
			if((lbm.flags[n]&TYPE_S)==TYPE_S) continue;
			uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
			const uchar fl = lbm.flags[n];
			const double ux=(double)lbm.u.x[n], uy=(double)lbm.u.y[n], uz=(double)lbm.u.z[n];
			const double umag = sqrt(sq(ux)+sq(uy)+sq(uz));
			if(umag>max_u) { max_u=umag; mux=x; muy=y; muz=z; muf=fl; }
			const double drho = fabs((double)lbm.rho[n]-1.0);
			if(drho>max_drho) max_drho = drho;
			if(lbm.T[n]>max_T_all) { max_T_all=lbm.T[n]; mTx=x; mTy=y; mTz=z; mTf2=fl; }
			const bool ordinary = (fl&(TYPE_E|TYPE_T))==0u;
			if(ordinary) { min_Tf=fmin(min_Tf, lbm.T[n]); max_Tf=fmax(max_Tf, lbm.T[n]); nTf++; }
			else         { min_Tb=fmin(min_Tb, lbm.T[n]); max_Tb=fmax(max_Tb, lbm.T[n]); nTb++; }
			const bool in_slot = fabs((float)y+0.5f-yc)<=lbm_inlet_half_w;
			if(z==Nz-1u && in_slot && (fl&TYPE_E)) sum_in_u += -uz;
			if(z==Nz-3u && in_slot && ordinary) { sum_in_p += (double)units.p_from_rho(lbm.rho[n]); sum_near_T += (double)lbm.T[n]; sum_in_uz_real += -uz; p_in_n++; near_n++; }
			if(y==1u || y==Ny-2u) {
				sum_out_p += (double)units.p_from_rho(lbm.rho[n]); p_out_n++;
				const double outward = y==1u ? -uy : uy;
				if(outward>0.0) { sum_out_uT += outward*(double)lbm.T[n]; sum_out_u += outward; }
			}
		}
		const bool blown_up = max_u!=max_u || max_drho!=max_drho || max_drho>0.45;
		const double Tout_lbm = sum_out_u>0.0 ? sum_out_uT/sum_out_u : (double)lbm_T_in;
		const float si_T_out  = units.si_T((float)Tout_lbm);
		const float si_T_near = near_n ? units.si_T((float)(sum_near_T/(double)near_n)) : si_T_in;
		const float si_dT     = si_T_out-si_T_in;
		const float si_P_strip = si_rho*(units.si_u((float)sum_out_u)*si_dx*si_dx)*si_cp*si_dT; // heat removed by MEASURED outflow (honest while flow develops)
		const float si_P_full  = si_P_strip*si_fin_len_x/si_strip_w;
		const float si_dp = (p_in_n&&p_out_n) ? units.si_p((float)(sum_in_p/(double)p_in_n))-units.si_p((float)(sum_out_p/(double)p_out_n)) : 0.0f;
		const double balance = sum_in_u>0.0 ? sum_out_u/sum_in_u : 0.0;

		// adaptive inlet pressure: a TYPE_E inlet imposes rho AND u; with rho=1 at both inlet and outlets the net driving
		// pressure is zero and the flow stalls below the target rate (mass out/in << 1). Raise the preset inlet density
		// each block until the measured outflow matches the target. Steady value = the true pump head for this geometry.
		drho_in = fmax(0.0, fmin(drho_in+0.03*(1.0-balance), 0.15)); // 3x gain: reaches the operating point in ~30-50k steps
		parallel_for(lbm.get_N(), [&](ulong n) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
			if(z==Nz-1u && (lbm.flags[n]&(TYPE_E|TYPE_T))==(TYPE_E|TYPE_T)) lbm.rho[n] = 1.0f+(float)drho_in;
		});
		lbm.rho.write_to_device();

		print_info("t="+to_string(lbm.get_t())
			+" | U_in="+to_string(units.si_u((float)(near_n ? sum_in_uz_real/(double)near_n : 0.0)), 3)+"/"+to_string(si_u_in, 3)+" m/s"
			+" | p_in="+to_string(units.si_p((float)(drho_in/3.0))/1000.0f, 3)+" kPa"
			+" | U_max="+to_string(units.si_u((float)max_u), 3)+" m/s @("+to_string(mux)+","+to_string(muy)+","+to_string(muz)+") f="+to_string((uint)muf)
			+" | mass out/in="+to_string((float)balance, 3)
			+" | max|rho-1|="+to_string((float)max_drho, 4));
		print_info("   dP="+to_string(si_dp/1000.0f, 3)+" kPa"
			+" | T_near_in="+to_string(si_T_near-273.15f, 2)
			+" | T_out="+to_string(si_T_out-273.15f, 2)+" degC (dT="+to_string(si_dT, 2)+" K)"
			+" | Q_strip="+to_string(si_P_strip, 3)+" W | Q_array~"+to_string(si_P_full, 1)+" W (upper bound)");
		print_info("   T fluid=["+to_string(units.si_T(min_Tf)-273.15f, 1)+","+to_string(units.si_T(max_Tf)-273.15f, 1)+"] degC n="+to_string(nTf)
			+" | T bound=["+to_string(units.si_T(min_Tb)-273.15f, 1)+","+to_string(units.si_T(max_Tb)-273.15f, 1)+"] n="+to_string(nTb)
			+" | T_max @("+to_string(mTx)+","+to_string(mTy)+","+to_string(mTz)+") f="+to_string((uint)mTf2));
		if(units.si_T(max_Tf)>si_T_wall+2.0f) print_warning("thermal overshoot above T_wall: D3Q7 dispersion (tau_T near 0.5); results near walls less reliable; refine cells_per_pitch if it grows.");
		if(blown_up) { print_warning("SIMULATION UNSTABLE (NaN or |rho-1|>0.45). Lower lbm_u (raises pressure unit) and restart. Stopping."); wait(); break; }

		if(export_vtk && lbm.get_t()>=next_vtk) {
			lbm.rho.write_device_to_vtk(vtk_path); lbm.u.write_device_to_vtk(vtk_path);
			lbm.T.write_device_to_vtk(vtk_path); lbm.flags.write_device_to_vtk(vtk_path, false);
			next_vtk += vtk_every;
		}

		conv_samples++;
		if(conv_samples>=3u && lbm.get_t()>=t_min && balance>0.95 && converged(Tout2, Tout1, (double)si_T_out, 1E-5)) {
			print_info("Converged on outlet temperature.");
			print_info("RESULT: dP = "+to_string(si_dp/1000.0f, 3)+" kPa | T_out = "+to_string(si_T_out-273.15f, 2)+" degC (dT = "+to_string(si_dT, 2)+" K) | Q_array <= "+to_string(si_P_full, 1)+" W at isothermal "+to_string(si_T_wall-273.15f, 0)+" degC wall");
			wait();
			break;
		}
		Tout2 = Tout1; Tout1 = (double)si_T_out;
	}
} /**/

defines.hpp

//#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 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 TEMPERATURE_LIMIT 0.11f // optional D3Q7 dispersion limiter: clamps temperature to T_avg +- this value (LBM units) every step; set slightly above the physical T range mapped via units.set_m_kg_s_K(); comment out to disable (required for setups like Rayleigh-Benard whose T range exceeds it)
//#define SUBGRID // disabled: channel Re ~ 90 (laminar) in the fin-array setup; 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.0001f // pressure/density view sensitivity; rho color range is [1.0f-GRAPHICS_RHO_DELTA, 1.0f+GRAPHICS_RHO_DELTA]
#define GRAPHICS_T_DELTA 0.10f // 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 4u // set how many streamlines there are every x lattice points (lowered 8->4 for denser streamlines in the narrow channels)
#define GRAPHICS_STREAMLINE_LENGTH 128u // 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_LBS 8u // local box size for local memory optimization in graphics_q() kernel, possible values: 0u (disable local memory optimization), 4u (no speedup), 8u (default, ~10% speedup)

//#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