Cells excluded from the simulation box ?

[O-Schell]

First of all: as someone just starting in CFD, I’m continually amazed by how seamlessly your code balances power and clarity. Truly, I can’t thank you enough for sharing this with the world !

I'm fully simulating a cyclorotor but when activating MOVING_BOUNDARIES, it seems to remove from the box of simulation the cells that were part of the STL.

I've try counted the cells that no flag but the number is steady; I tried settings all the cells to G after unvoxeling the stl but didn't change a thing; I've refactored my code in case my badly styled code was hiding a mistake and spend many hours testing more or less anything coming to my mind that could potentially influence what I was seeing. And yet, at my despair, I cannot tell you exactly is going wrong

Some information about my code
I have FORCE_FIELDS (although removed for the screen), EQUILIBRIUM_BOUNDARIES, MOVING_BOUNDARIES and SUBGRID activated. Like in most of the examples (which are amazing to learn and start), for the moving boundaries, I'm setting the instantaneous speed in the voxelization.

Does anyone have an idea of where the problem might be ? I haven't found anything that could explain this behavior and any lead would help a long way since I'm running of where to look / system to debug

Here's my code and below an illustration of the problem:

// ##################################################################
// Full cinematic for 4 Blades; required extensions in defines.hpp: FP16S, EQUILIBRIUM_BOUNDARIES, MOVING_BOUNDARIES, SUBGRID, INTERACTIVE_GRAPHICS or GRAPHICS
// ##################################################################
void main_setup() {
#include <fstream>
#include <iostream>
#include <cmath>
#include <string>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

	// ##################################################################
	// Setup: CSV, parameters, and units
	// ##################################################################
	std::ofstream csv_file("simulation_data__full_cinematic.csv");
	if (!csv_file.is_open()) {
		print_error("Failed to open CSV file for writing!");
		return;
	}
	csv_file << "t,flap_angle_deg,force_x,force_y,force_z,torque_x,torque_y,torque_z\n";

	// --- Blade parameters ---
	const float blade_A_radius = 0.10f;          // [m] Distance to center of rotation A
	const float blade_span = 0.15f;              // [m] Blade span
	const float blade_chord = 0.05f;             // [m] Blade chord
	const float blade_thickness = 0.15f * blade_chord; // [m] Thickness (15% of chord)
	const float blade_A_rotation_speed = 1100;   // [RPM] Rotational velocity
	const float blade_B_max_angle = 20.0f;      // [deg] Max flap angle
	const float blade_B_max_angle_rad = blade_B_max_angle * M_PI / 180.0f;
	const float margin_box_size = 0.8f;
	const float box_vertical_size = 2 * margin_box_size * (blade_A_radius + blade_chord);
	const float box_horizontal_ratio = 1.1f * (blade_chord / (box_vertical_size / 2));
	const int number_of_turns = 5;

	// --- Simulation units and LBM setup ---
	const uint3 lbm_N = resolution(float3(box_horizontal_ratio, 1.5f, 1.5f), 500u);
	const float si_u = (2 * M_PI * blade_A_radius) * (blade_A_rotation_speed / 60.0f); // [m/s]
	const float si_length = blade_chord;         // [m] Characteristic length
	const float si_nu = 1.821e-5f;              // [m²/s] Kinematic viscosity (air @20°C)
	const float si_rho = 1.2041f;               // [kg/m³] Density (air @20°C)
	const float si_t = 1.0f / (blade_A_rotation_speed / 60.0f); // [s] Time for 1 rotation
	const float lbm_u = 0.075f;                 // LBM ideal velocity
	const float lbm_rho = 1.0f;                 // LBM ideal density

	units.set_m_kg_s(lbm_N.x, lbm_u, lbm_rho, box_horizontal_ratio * box_vertical_size, si_u, si_rho);

	const float lbm_length = units.x(blade_A_radius);
	const float lbm_omega = lbm_u / lbm_length;
	const float lbm_L = units.x(si_length);
	const float lbm_v = units.nu(si_nu);
	const float lbm_Re = (lbm_u * lbm_L) / lbm_v;
	const ulong lbm_T = units.t(si_t);
	const ulong lbm_dt = 1ull;

	LBM lbm(lbm_N, 1u, 1u, 1u, lbm_v);
	print_info(to_string(units.si_Re(lbm_Re)) + " = ( " + to_string(units.si_u(lbm_u)) + " * " + to_string(units.si_x(lbm_L)) + " ) / " + to_string(units.si_nu(lbm_v)));
	print_info("Re = " + to_string(to_uint(units.si_Re(si_length, si_u, si_nu))));
	print_info(to_string(si_t, 3u) + " seconds = " + to_string(lbm_T) + " time steps");

	// ##################################################################
	// Define geometry: Load 4 blades (0°, 90°, 180°, 270°)
	// ##################################################################
	float3x3 init_rot(float3(0, 0, 1), 3 * M_PI / 2);
	Mesh* meshes[4];
	float3 centers[4] = {
		float3(lbm.center().x, lbm.center().y, lbm.center().z),
		float3(lbm.center().x, lbm.center().y, lbm.center().z),
		float3(lbm.center().x, lbm.center().y, lbm.center().z),
		float3(lbm.center().x, lbm.center().y, lbm.center().z)
	};

	for (int i = 0; i < 4; ++i) {
		meshes[i] = read_stl(get_exe_path() + "../stl/NACA0015_L50_W150.stl", lbm.size(), centers[i], init_rot, lbm_length);
		if (!meshes[i]) {
			print_error("Failed to load mesh for blade " + to_string(i));
			return;
		}
	}

	// ##################################################################
	// Boundary conditions
	// ##################################################################
	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 (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) {
			lbm.flags[n] = TYPE_E;
			lbm.rho[n] = 1.0f;
			lbm.u.x[n] = lbm.u.y[n] = lbm.u.z[n] = 0.0f;
		}
		if (lbm.flags[n] != TYPE_S) {
			lbm.u.y[n] = units.u(1.0f / 3.6f); // km/h → m/s → lattice units
		}
		});

	// ##################################################################
	// Initialize simulation
	// ##################################################################
	lbm.run(0u, number_of_turns * lbm_T);

	print_info("BOX SIZE (LAT): " + to_string(lbm_N.x) + ", " + to_string(lbm_N.y) + ", " + to_string(lbm_N.z));
	print_info("SIM TIME LENGTH (LAT): " + to_string(number_of_turns * si_t, 4) + "s (" + to_string(number_of_turns * lbm_T) + " lattice units)");

	// --- Blade state tracking ---
	float rot_angles[4] = { 0 };
	float trans_x[4] = { 0 }, trans_y[4] = { 0 };

	// ##################################################################
	// Main simulation loop
	// ##################################################################
	while (lbm.get_t() < number_of_turns * lbm_T) {
		const float t = lbm.get_t();
		const float current_si_t = units.si_t(t);
		const float rad_rotation_progress = float_mod(2 * M_PI * current_si_t / si_t, 2 * M_PI);

		// --- Update each blade ---
		for (int i = 0; i < 4; ++i) {
			const float phase_offset = i * M_PI / 2.0f;
			const float rad_progress = float_mod(rad_rotation_progress + phase_offset, 2 * M_PI);

			// Flap angle: sync + oscillation
			const float flap_angle_sync = float_mod(rad_progress + (M_PI / 2), 2 * M_PI);
			const float flap_angle_osci = -sin(rad_progress) * blade_B_max_angle_rad;
			const float flap_angle = float_mod(flap_angle_sync + flap_angle_osci, 2 * M_PI);
			const float rotation_needed = flap_angle - rot_angles[i];

			// Translation
			const float Bx = blade_A_radius * cos(rad_progress);
			const float By = blade_A_radius * sin(rad_progress);
			const float trans_needed_x = Bx - trans_x[i];
			const float trans_needed_y = By - trans_y[i];

			// Apply transformations
			lbm.unvoxelize_mesh_on_device(meshes[i], TYPE_S);
			meshes[i]->rotate(float3x3(float3(1.0f, 0.0f, 0.0f), rotation_needed));
			meshes[i]->translate(float3(0.0f, units.x(trans_needed_x), units.x(trans_needed_y)));
			lbm.voxelize_mesh_on_device(meshes[i], TYPE_S, meshes[i]->get_center_of_mass(), float3(0.0f), float3(0.0f, lbm_omega, 0.0f));

			// Update state
			rot_angles[i] = flap_angle;
			trans_x[i] = Bx;
			trans_y[i] = By;
		}

		// --- Log forces/torques to CSV ---
		//const float3 lbm_com = lbm.object_center_of_mass(TYPE_S);
		//const float3 lbm_force = lbm.object_force(TYPE_S);
		//const float3 lbm_torque = lbm.object_torque(lbm_com, TYPE_S);
		//csv_file << current_si_t << ","
		//	<< (rad_rotation_progress * 180.0f / M_PI) << ","
		//	<< units.si_F(lbm_force.x) << "," << units.si_F(lbm_force.y) << "," << units.si_F(lbm_force.z) << ","
		//	<< units.si_M(lbm_torque.x) << "," << units.si_M(lbm_torque.y) << "," << units.si_M(lbm_torque.z) << "\n";

		lbm.run(lbm_dt, lbm_T);

#if defined(GRAPHICS) && !defined(INTERACTIVE_GRAPHICS)
		if (lbm.graphics.next_frame(lbm_T, 30.0f)) {
			lbm.graphics.set_camera_free(
				float3(0.353512f * (float)Nx, -0.150326f * (float)Ny, 1.643939f * (float)Nz),
				-25.0f, 61.0f, 100.0f
			);
			lbm.graphics.write_frame();
		}
#endif
	}
	csv_file.close();
}

[O-Schell]

After more time, it seems that the main culprit is the rotation_center paramter in voxelize_mesh_on_device; It seems that when the values are too far off the center of the simulation box (but still well inside it), it produce this bug. I'm still investigating

[O-Schell]

Ok, so as far as I understand (with the help of Claude), the voxelize_mesh_on_device doesn't support a moving rotation_center (the nominal use case is pure rotating object): It's using a hash of the cells to do some calculation but the problem is that if the rotation center move, the hash has changed and unexpected behavior arise...
The concerned code is in kernel.cpp around line 2337 (version 3.7)

[ProjectPhysX]

Hi @O-Schell,

voxelization applies an instantaneous velocity to all voxelized cells, that can be from rotation or translation or a combination of both. There is no hash of cells. Between voxelization steps you can also modify/move the rotation center.
Please don't rely on LLM hallucinations.
Will look into your issue in detail some time, please be patient.

Kind regards,
Moritz