Galilean Invariance and Force Calculation on Moving Boundaries in FluidX3D

Thanks for your incredible work on FluidX3D! As a beginner in this field, I am deeply impressed by the computation speed and the clean architecture of this code. It has provided me with an exceptional learning experience and a powerful tool for my research.
Recently, I have been conducting some benchmarks using moving meshes (specifically a NACA airfoil case) to validate Galilean Invariance within the framework. I have attached my main_setup.cpp for your reference. In the code, I implemented a mode switch at the beginning of the main_setup function to toggle between two physically equivalent setups:

// 1 = Wind Tunnel Mode (Stationary mesh, moving fluid)
// 2 = Free Flight Mode (Moving mesh, stationary fluid)
const int TEST_MODE = 1;

I have been logging force results and plotted a comparison of the lift and drag coefficients for both modes. However, the results revealed some discrepancies that I am struggling to reconcile:

1. Lift Coefficient: While the results are positive in both modes, the magnitudes are not equal. This suggests a potential discrepancy in the force integration or the way moving boundaries are handled during the voxelization/streaming steps.
2. Drag Coefficient: Most notably, in "Free Flight Mode," the calculated drag is negative (effectively acting as a thrust). This is physically incorrect for a passive gliding case and points to a lack of Galilean invariance in the current force calculation.

I was wondering if you might have any insights into this phenomenon regarding Galilean Invariance with moving boundaries in FluidX3D, or if this is a known challenge that others have encountered. Specifically:

1. How does FluidX3D currently handle aerodynamic force calculations on moving boundaries?
2. If the current implementation is optimized primarily for static meshes, could you suggest which theoretical framework or literature I should refer to for implementing a correction?
3. In which kernel functions would you recommend I implement these modifications to best align with the FluidX3D architecture?
   Thank you very much for your time and for sharing this powerful tool with the community. I look forward to your professional insights.

#include "setup.hpp"
#include <fstream>
#include <iostream>
#include <iomanip>
#include <limits>
#include <cmath>
#include <filesystem>
#include "units.hpp"

//Comparison of lift and drag coefficients between stationary airfoil with inflow and moving airfoil without inflow

// --------------------------------------------------------
// Geometry Generation Tool (NACA Airfoil)
// --------------------------------------------------------
Mesh* create_naca_mesh(float chord, float thickness, float z_width, float3 center_offset) {
	const int N_points = 100;
	std::vector<float3> profile_top(N_points), profile_bot(N_points);

	// generate airfoil coordinate points
	for (int i = 0; i < N_points; ++i) {
		float pos = (float)i / (float)(N_points - 1);
		// NACA 4-digit airfoil thickness formula
		float yt = 5.0f * thickness * (0.2969f * sqrt(pos) - 0.1260f * pos - 0.3516f * pos * pos + 0.2843f * pos * pos * pos - 0.1015f * pos * pos * pos * pos);
		float x = (pos - 0.5f) * chord;
		float y = yt * chord;
		profile_top[i] = float3(x, y, 0.0f);
		profile_bot[i] = float3(x, -y, 0.0f);
	}
	// force trailing edge closure (prevent floating point errors)
	profile_top[N_points - 1].y = 0.0f;
	profile_bot[N_points - 1].y = 0.0f;

	float z0 = -0.5f * z_width, z1 = +0.5f * z_width;

	Mesh* mesh = new Mesh((N_points - 1) * 8, center_offset);
	int t_idx = 0;

	// add a quad (composed of two triangles)
	auto add_quad = [&](float3 p1, float3 p2, float3 p3, float3 p4) {
		mesh->p0[t_idx] = p1; mesh->p1[t_idx] = p2; mesh->p2[t_idx] = p3; t_idx++;
		mesh->p0[t_idx] = p3; mesh->p1[t_idx] = p2; mesh->p2[t_idx] = p4; t_idx++;
	};

	for (int i = 0; i < N_points - 1; ++i) {
		// top surface
		add_quad(
			float3(profile_top[i].x, profile_top[i].y, z0),
			float3(profile_top[i + 1].x, profile_top[i + 1].y, z0),
			float3(profile_top[i].x, profile_top[i].y, z1),
			float3(profile_top[i + 1].x, profile_top[i + 1].y, z1)
		);

		// bottom surface
		add_quad(
			float3(profile_bot[i + 1].x, profile_bot[i + 1].y, z0),
			float3(profile_bot[i].x, profile_bot[i].y, z0),
			float3(profile_bot[i + 1].x, profile_bot[i + 1].y, z1),
			float3(profile_bot[i].x, profile_bot[i].y, z1)
		);

		// left cap (Z = z0 plane) - close the side
		add_quad(
			float3(profile_top[i].x, profile_top[i].y, z0),
			float3(profile_top[i + 1].x, profile_top[i + 1].y, z0),
			float3(profile_bot[i].x, profile_bot[i].y, z0),
			float3(profile_bot[i + 1].x, profile_bot[i + 1].y, z0)
		);

		//  right cap (Z = z1 plane) - close the side
		add_quad(
			float3(profile_top[i].x, profile_top[i].y, z1),
			float3(profile_bot[i].x, profile_bot[i].y, z1),
			float3(profile_top[i + 1].x, profile_top[i + 1].y, z1),
			float3(profile_bot[i + 1].x, profile_bot[i + 1].y, z1)
		);
	}

	mesh->find_bounds(); 
	return mesh;
}

// Define Motion State Structure
struct MotionState {
	float3 velocity;    // linear velocity (for fluid calculation)
	float3 omega;     // angular velocity (for fluid calculation)
	float3 d_pos;       // displacement for the current frame (for mesh movement)
	float  d_theta;     // rotation increment for the current frame (radians, for mesh rotation)
};


// translation
// to simulate free flight: stationary fluid, airfoil flying at constant speed
MotionState calc_pure_translation(float speed_x) {
	MotionState state;

	state.velocity = float3(-speed_x, 0.0f, 0.0f);
	state.omega = float3(0.0f);
	state.d_pos = state.velocity; // dt=1
	state.d_theta = 0.0f;
	return state;
}


void main_setup() { // required extensions in defines.hpp: D2Q9, FP16S, MOVING_BOUNDARIES, SUBGRID, FORCE_FIELD, INTERACTIVE_GRAPHICS or GRAPHICS

	// ===================================== [Mode Switch] =====================================
	// 1 = Wind Tunnel Mode (Stationary mesh, moving fluid) -> Traditional CFD
	// 2 = Free Flight Mode (Moving mesh, stationary fluid) -> Galilean Invariance validation
	// =========================================================================================
	const int TEST_MODE = 1;

	// parameter settings
	const uint R = 16u;
	const float chord_lbm = 12.0f * (float)R; // characteristic length: chord L
	const float flow_speed_mag = 0.1f;        

	// determine inflow velocity (u_inf) and airfoil velocity (mesh_vel) based on mode
	float u_inf_lbm_set = 0.0f;    // velocity passed to LBM boundaries and initial field
	float mesh_speed_set = 0.0f;   // movement speed passed to the Mesh
	if (TEST_MODE == 1) {          // wind tunnel
		u_inf_lbm_set = flow_speed_mag;
		mesh_speed_set = 0.0f;
	}
	else {                         // free flight: ambient air is stationary
		u_inf_lbm_set = 0.0f;         
		mesh_speed_set = flow_speed_mag; 
	}

	const float rho_lbm = 1.0f; // LBM density (rho)
	const float q_lbm = 0.5f * 1.0f * flow_speed_mag * flow_speed_mag; // LBM dynamic pressure

	// define physical (SI) parameters
	const float chord_si = 0.1f;                            // chord length [m]
	const float u_inf_si = 0.01f;                           // flow velocity [m/s]
	const float rho_si = 1.225f;                            // density: air 1.225 [kg/m^3]
	const float mu_si = 1.79e-5f;                           // dynamic viscosity: air 1.79e-5 [Pa s]
	const float nu_si = mu_si / rho_si;                     // kinematic viscosity [m^2/s]
	const float Re_si = u_inf_si * chord_si / nu_si;        // Reynolds number; LBM Reynolds number is set equal to this
	const float q_si = 0.5f * rho_si * u_inf_si * u_inf_si; // dynamic pressure

	Units units; 

	// mapping
	units.set_m_kg_s(
		chord_lbm, flow_speed_mag, rho_lbm,  
		chord_si, u_inf_si, rho_si    
	);

	// initialize LBM
	const uint Nx = 256u * R;
	const uint Ny = 128u * R;
	const uint Nz = 1u;
	LBM lbm(Nx, Ny, Nz, units.nu_from_Re(Re_si, chord_lbm, flow_speed_mag)); 

	// geometry initialization
	Mesh* mesh = create_naca_mesh(chord_lbm, 0.12f, 10.0f, float3(0.0f));

	// 1. initial Attitude
	float initial_angle_deg = -5.0f;
	float initial_angle_rad = radians(initial_angle_deg);
	mesh->rotate(float3x3(float3(0, 0, 1), initial_angle_rad));

	// adjust starting position based on mode
	float3 start_pos;
	if (TEST_MODE == 1) {
		// wind tunnel mode: place on the left, fluid blows from left to right
		start_pos = float3(Nx * 0.3f, Ny * 0.5f, 0.5f);
	}
	else {
		// free flight: place on the right, airfoil flies from right to left (to ensure sufficient "runway")
		start_pos = float3(Nx * 0.8f, Ny * 0.5f, 0.5f);
	}

	float3 current_pos = start_pos; 
	float current_angle_rad = initial_angle_rad; 
	mesh->translate(start_pos);     

	//lbm.voxelize_mesh_on_device(mesh, TYPE_S, mesh->get_center());

	// initial voxelization
	MotionState init_state = calc_pure_translation(mesh_speed_set);
	lbm.voxelize_mesh_on_device(mesh, TYPE_S, mesh->get_center(), init_state.velocity, init_state.omega);

	// boundary conditions
	parallel_for(lbm.get_N(), [&](ulong n) {
		uint x, y, z;
		lbm.coordinates(n, x, y, z);
		if (x == 0 || x == Nx - 1 || y == 0 || y == Ny - 1) { // boundary
			lbm.flags[n] = TYPE_E;
			lbm.u.x[n] = u_inf_lbm_set;
		}
		else {                                                // non-boundary
			if (lbm.flags[n] != TYPE_S) { 
				lbm.u.x[n] = u_inf_lbm_set; 
			}
		}
		});
	lbm.run(0u); 

	// initialize flow field and upload
	lbm.graphics.visualization_modes = VIS_FLAG_SURFACE | VIS_FIELD; 

	const ulong total_time = 1000000ull;

	// ===================== Output Settings =====================
	std::filesystem::create_directories(get_exe_path() + "export/");
	const std::string outfile = get_exe_path() + "export/forces.csv";
	std::ofstream fout(outfile);
	fout << "t_step,t_si,Fx_lbm,Fy_lbm,Mz_lbm,L_lbm,D_lbm,CL_lbm,CD_lbm,Cmz_lbm,Fx_si,Fy_si,Mz_si,L_si,D_si,CL_si,CD_si,Cmz_si\n";
	fout << std::setprecision(10);

	if (!fout.is_open()) {
		print_info("ERROR: cannot open " + outfile);
		return;
	}

	const uint  sample_every = 10u;      // sampling interval

	while (lbm.get_t() < total_time)
	{
		ulong t = lbm.get_t();

		MotionState state;
		if (TEST_MODE == 1) {   // wind tunnel: static mesh
			state = calc_pure_translation(0.0f);
		}
		else {                  // free flight: mesh in uniform motion
			state = calc_pure_translation(mesh_speed_set);
			if (current_pos.x < chord_lbm * 2.0f) {
				print_info("Airfoil reached left boundary. Stopping.");
				break;
			}
		}

		lbm.run(10u);
		state.d_pos = state.velocity * 10.f;
		state.d_theta = state.omega.z * 10.f;

		if (t % sample_every == 0u) {
			lbm.update_force_field();
			const float3 F_lbm = lbm.object_force(TYPE_S);

			const float Fx_lbm = F_lbm.x;
			const float Fy_lbm = F_lbm.y;
			const float L_lbm = Fy_lbm;
			const float D_lbm = Fx_lbm;
			const float CL_lbm = L_lbm / (q_lbm * chord_lbm * 1.0f); // use unit span for 2D airfoil
			const float CD_lbm = D_lbm / (q_lbm * chord_lbm * 1.0f); // use unit span for 2D airfoil

			// calculate offset vector from 1/4 chord point relative to the mesh center
			// mesh center is at 0.5c, target is at 0.25c, so offset "forward" by 0.25c
			float offset_dist = -0.25f * chord_lbm;

			// rotate the offset vector (around Z-axis)
			float3 offset_vec;
			offset_vec.x = offset_dist * cos(current_angle_rad);
			offset_vec.y = offset_dist * sin(current_angle_rad);
			offset_vec.z = 0.0f;

			// obtain the absolute coordinates of the 1/4 chord point
			float3 center_quarter_chord = mesh->get_center() + offset_vec;

			// calculate moment [LBM units]
			float3 Torque_lbm_vec = lbm.object_torque(center_quarter_chord, TYPE_S);
			const float Mz_lbm = Torque_lbm_vec.z;

			// calculate moment coefficient Cm = Moment / (q * S * c)
			const float Cmz_lbm = Mz_lbm / (q_lbm * chord_lbm * 1.0f * chord_lbm);

			// LBM -> SI
			const float t_si = units.si_t((float)lbm.get_t()); // [s] 
			const float Fx_si = units.si_F(Fx_lbm);   // [N]
			const float Fy_si = units.si_F(Fy_lbm);   // [N]
			const float Mz_si = units.si_M(Mz_lbm);   // [Nm]

			// Lift/Drag for freestream along +X:
			const float L_si = Fy_si;   // [N/m]
			const float D_si = Fx_si;   // [N/m]
			const float CL_si = L_si / (q_si * chord_si * units.si_x(1.0f)); 
			const float CD_si = D_si / (q_si * chord_si * units.si_x(1.0f)); 
			const float Cmz_si = Mz_si / (q_si * chord_si * units.si_x(1.0f) * chord_si); // M / (q * S * chord)

			fout << lbm.get_t() << "," << t_si << ","
				<< Fx_lbm << "," << Fy_lbm << "," << Mz_lbm << ","
				<< L_lbm << "," << D_lbm << "," << CL_lbm << "," << CD_lbm << "," << Cmz_lbm << ","
				<< Fx_si << "," << Fy_si << "," << Mz_si << ","
				<< L_si << "," << D_si << "," << CL_si << "," << CD_si << "," << Cmz_si << "\n";
		}
	
		lbm.unvoxelize_mesh_on_device(mesh, TYPE_F);

		mesh->rotate(float3x3(float3(0, 0, 1), state.d_theta));
		mesh->translate(state.d_pos);

		lbm.voxelize_mesh_on_device(mesh, TYPE_S, mesh->get_center(), state.velocity, state.omega);

		current_angle_rad += state.d_theta;
		current_pos += state.d_pos;
	}
	fout.flush();
	fout.close();
}

Start with the attached main_setup.cpp and reproduce both TEST_MODE values, recording the generated forces.csv output. Then trace voxelize_mesh_on_device, update_force_field, object_force, and object_torque, along with setup.hpp and units.hpp, to understand the moving-boundary force path. Done means explaining the discrepancy and identifying a validated correction or clearly documenting the limitation.