ALM for VAWT Simulation Issues and GPU Optimization in FluidX3D

[IsTrivial]

Hello everyone. I’ve spent the last few days trying to simulate a 2D hydrokinetic VAWT using the Actuator Line Model (ALM) and FluidX3D’s FORCE_FIELD, driven by my own aerodynamic coefficient data (XFOIL). The physics implementation is going well (L/D calculation, dynamic stall, etc.), but I’m hitting a huge bottleneck because I have to copy the entire fluid velocity field (lbm.u) from the GPU to the CPU at every step in order to perform interpolation and compute the force on the blade. I was wondering if there is any internal function or utility that would let me sample only the interpolated velocity right at the actuator-line location and keep most of the computation on the GPU. Also, looking ahead to moving this to 3D, keeping the same cell resolution across the whole domain will be unfeasible. Is there any supported or in-development way to implement an adaptive mesh so that high resolution can be concentrated only around the rotor? Any advice on how other users have approached integrating lookup tables into the FluidX3D kernel would be very helpful.

#include "setup.hpp"
#include <fstream>
#include <vector>
#include <deque>
#include <stdexcept>
#include <cmath>
#include <iostream>
#include <iomanip>

using namespace std;

// ============================================================================
// ESTRUCTURAS
// ============================================================================
struct BladeState {
    float alpha_filtered = 0.0f;
    float alpha_prev = 0.0f;
    float Un_prev = 0.0f;
    deque<float> da_dt_buffer;
};

// ============================================================================
// INTERPOLACIÓN Y UTILIDADES
// ============================================================================
float3 get_interpolated_velocity_cpu(LBM& lbm, float3 pos) {
    int ix = clamp((int)floor(pos.x), 0, (int)lbm.get_Nx() - 2);
    int iy = clamp((int)floor(pos.y), 0, (int)lbm.get_Ny() - 2);

    float fx = pos.x - ix;
    float fy = pos.y - iy;

    auto get_vec = [&](int x, int y) {
        ulong n = lbm.index(x, y, 0);
        return float3(lbm.u.x[n], lbm.u.y[n], lbm.u.z[n]);
        };

    float3 u00 = get_vec(ix, iy);     float3 u01 = get_vec(ix, iy + 1);
    float3 u10 = get_vec(ix + 1, iy); float3 u11 = get_vec(ix + 1, iy + 1);

    return (u00 * (1.0f - fy) + u01 * fy) * (1.0f - fx) + (u10 * (1.0f - fy) + u11 * fy) * fx;
}

float3 sample_velocity_lineAverage(LBM& lbm, float3 center, float radius, float3 unit_x, float3 unit_y, int N_points = 80) {
    float3 v_avg(0.0f);
    int valid = 0;
    for (int i = 0; i < N_points; i++) {
        float theta = 2.0f * pi * i / N_points; 
        float3 sample_pos = center + radius * (cos(theta) * unit_x + sin(theta) * unit_y);

        if (sample_pos.x >= 0 && sample_pos.x < lbm.get_Nx() && sample_pos.y >= 0 && sample_pos.y < lbm.get_Ny()) {
            v_avg += get_interpolated_velocity_cpu(lbm, sample_pos);
            valid++;
        }
    }
    return valid > 0 ? v_avg / (float)valid : float3(0.0f);
}

void smear_force(LBM& lbm, const float3& pos, const float3& f_fluid, float epsilon_lu) {
    int range = (int)ceilf(3.0f * epsilon_lu);
    float sigma_sq = epsilon_lu * epsilon_lu;
    float gauss_norm = 1.0f / (2.0f * pi * sigma_sq); 
    int zz = 0; 

    int min_x = max(0, (int)floor(pos.x) - range);
    int max_x = min((int)lbm.get_Nx() - 1, (int)floor(pos.x) + range);
    int min_y = max(0, (int)floor(pos.y) - range);
    int max_y = min((int)lbm.get_Ny() - 1, (int)floor(pos.y) + range);

    for (int xx = min_x; xx <= max_x; ++xx) {
        for (int yy = min_y; yy <= max_y; ++yy) {
            float3 dr = float3(xx + 0.5f, yy + 0.5f, zz + 0.5f) - pos;
            float r2 = dot(dr, dr);

            if (r2 > 9.0f * sigma_sq) continue;

            float g = gauss_norm * expf(-r2 / (2.0f * sigma_sq));
            ulong n = lbm.index(xx, yy, zz);

            lbm.F.x[n] += f_fluid.x * g;
            lbm.F.y[n] += f_fluid.y * g;
        }
    }
}

// ============================================================================
// CARGA DE DATOS 
// ============================================================================
bool read_airfoil_data(const string& filename, vector<float>& alphas, vector<float>& Re_nums,
    vector<vector<float>>& cls, vector<vector<float>>& cds, vector<vector<float>>& cms) {
    ifstream file(filename);
    if (!file) return false;
    string line, token;
    if (!getline(file, line)) return false;
    stringstream ss(line);
    getline(ss, token, ',');
    while (getline(ss, token, ',')) {
        token.erase(token.find_last_not_of(" \n\r\t") + 1);
        if (token == "|" || token.empty()) break;
        alphas.push_back(stof(token));
    }
    while (getline(file, line)) {
        stringstream row(line);
        getline(row, token, ',');
        if (token.empty()) continue;
        try { Re_nums.push_back(stof(token)); }
        catch (...) { continue; }
        vector<float> cl_row, cd_row, cm_row;
        int section = 0;
        while (getline(row, token, ',')) {
            token.erase(token.find_last_not_of(" \n\r\t") + 1);
            if (token == "|") { section++; continue; }
            if (token.empty()) continue;
            float val = stof(token);
            if (section == 0) cl_row.push_back(val);
            else if (section == 1) cd_row.push_back(val);
            else if (section == 2) cm_row.push_back(val);
        }
        cls.push_back(cl_row); cds.push_back(cd_row);
        cms.push_back(cm_row.empty() ? vector<float>(cl_row.size(), 0.0f) : cm_row);
    }
    return !Re_nums.empty();
}

float interpolate_2d(const vector<float>& alphas, const vector<float>& Re_nums,
    const vector<vector<float>>& values, float alpha, float Re) {
    if (alphas.empty()) return 0.0f;
    alpha = clamp(alpha, alphas.front(), alphas.back()); 
    Re = clamp(Re, Re_nums.front(), Re_nums.back());

    auto ait = lower_bound(alphas.begin(), alphas.end(), alpha);
    size_t ai = min((ait == alphas.begin()) ? 0 : (size_t)(ait - alphas.begin() - 1), alphas.size() - 2);
    float afrac = (alpha - alphas[ai]) / (alphas[ai + 1] - alphas[ai]);

    auto rit = lower_bound(Re_nums.begin(), Re_nums.end(), Re);
    size_t ri = min((rit == Re_nums.begin()) ? 0 : (size_t)(rit - Re_nums.begin() - 1), Re_nums.size() - 2);
    float Re_frac = (logf(Re) - logf(Re_nums[ri])) / (logf(Re_nums[ri + 1]) - logf(Re_nums[ri]));

    float v0 = values[ri][ai] + afrac * (values[ri][ai + 1] - values[ri][ai]);
    float v1 = values[ri + 1][ai] + afrac * (values[ri + 1][ai + 1] - values[ri + 1][ai]);
    return v0 + Re_frac * (v1 - v0);
}

// ============================================================================
// MAIN SETUP
// ============================================================================
void main_setup() {
    cout << unitbuf; cerr << unitbuf;
    cout << "=== INICIALIZANDO PARAMETROS ===" << endl;

    // --- PARÁMETROS FÍSICOS ---
    const float D_si = 1.0f;
    const float H_si = 1.0f;
    const float chord_si = 0.14f;
    const float U_inf = 1.0f;
    const float rpm = 36.3f;
    const int N_blades = 3;

    const float chord_pivot_ratio = 0.5f;
    const float AC_ratio = 0.25f;
    const float d_AC_BSC_si = (chord_pivot_ratio - AC_ratio) * chord_si;

    const float si_nu = 1.0e-6f;
    const float si_rho = 998.0f;
    const float R_si = D_si / 2.0f;
    const float si_omega = 2.0f * pi * (rpm / 60.0f); 
    const float si_period = 2.0f * pi / si_omega;

    // --- MALLA ---
    const float domain_x_si = 35.0f * D_si;
    const float domain_y_si = 20.0f * D_si;
    const uint3 lbm_N = resolution(float3(domain_x_si, domain_y_si, 0.0f), 1500u);

    // --- LBM SETUP ---
    const float lbm_u = 0.1f;
    units.set_m_kg_s((float)lbm_N.x, lbm_u, 1.0f, domain_x_si, U_inf, si_rho);

    const float lbm_nu = units.nu(si_nu);
    const float lbm_omega = units.omega(si_omega);
    const float lbm_R = units.x(R_si);
    const float chord_lbm = units.x(chord_si);
    const ulong lbm_T = units.t(5.0f * si_period);
    const ulong WARMUP_STEPS = 7000;


    LBM lbm(lbm_N, 1u, 1u, 1u, lbm_nu);

    // --- CONDICIONES DE FRONTERA ---
    parallel_for(lbm.get_N(), [&](ulong n) {
        lbm.rho[n] = 1.0f;
        lbm.u.x[n] = lbm_u; lbm.u.y[n] = 0.0f; lbm.u.z[n] = 0.0f;

        uint x, y, z; lbm.coordinates(n, x, y, z);
        if (x == 0u || x == lbm.get_Nx() - 1u || y == 0u || y == lbm.get_Ny() - 1u) {
            lbm.flags[n] = TYPE_E;
        }
        });

    lbm.rho.write_to_device();
    lbm.u.write_to_device();
    lbm.flags.write_to_device();

    // --- CARGA DE DATOS ---
    vector<float> alphas, Re_nums;
    vector<vector<float>> cls, cds, cms;
    if (!read_airfoil_data("airfoil_data.csv", alphas, Re_nums, cls, cds, cms)) return;

    // --- ARCHIVOS DE SALIDA ---
    ofstream torque_file("torque_data.txt");
    torque_file << "# Time\tTorque_Smooth" << endl;

    ofstream blade_file("blade_debug.csv");
    blade_file << "Time_s,Theta_deg,Phi_deg,Alpha_deg,Cl_3d,Cd_3d,Cn,Ct,V_mag,Torque_inst" << endl;

    // --- VARIABLES SIM ---
    vector<BladeState> blades(N_blades);
    const float alpha_smoothing_factor = 0.08f;
    const float3 center = float3(units.x(10.0f * D_si), lbm.center().y, 0.5f);

    float theta = 0.0f;
    deque<float> torque_output_buffer;
    double torque_sum = 0.0;
    long torque_samples = 0;
    const float wind_power_si = 0.5f * si_rho * (D_si * H_si) * pow(U_inf, 3.0f);

    cout << "=== STARTING MAIN LOOP ===" << endl;

    lbm.run(0u);

    while (lbm.get_t() < lbm_T) {
        parallel_for(lbm.get_N(), [&](ulong n) { lbm.F.x[n] = 0.0f; lbm.F.y[n] = 0.0f; });

        lbm.u.read_from_device();
        float torque_lbm_instant = 0.0f;

        for (int b = 0; b < N_blades; ++b) {
            float phi = theta + 2.0f * pi * b / N_blades;
            float phi_norm = fmod(phi, 2.0f * pi);
            if (phi_norm < 0) phi_norm += 2.0f * pi;

            float3 unit_r = float3(cos(phi), sin(phi), 0.0f);
            float3 unit_t = float3(-sin(phi), cos(phi), 0.0f);
            float3 pos_AC = center + lbm_R * unit_r;
            float3 blade_vel = lbm_omega * lbm_R * unit_t;

            float3 u_fluid = sample_velocity_lineAverage(lbm, pos_AC, 0.25f * chord_lbm, unit_t, unit_r, 60);
            float3 rel_vel = u_fluid - blade_vel;
            float V_mag = length(rel_vel); 

            if (V_mag < 0.001f) continue;

            float V_chord = dot(rel_vel, unit_t); 
            float V_normal = dot(rel_vel, unit_r);
            float alpha_rad = atan2(V_normal, -V_chord);

            // Corrección curvatura flujo
            alpha_rad += (si_omega * chord_si) / (4.0f * units.si_u(V_mag));
            float alpha_deg = degrees(alpha_rad); 

            float& alpha_f = blades[b].alpha_filtered;
            if (lbm.get_t() < 10) alpha_f = alpha_deg;
            alpha_f = alpha_f * (1.0f - alpha_smoothing_factor) + alpha_deg * alpha_smoothing_factor;

            float dt_si = units.si_t(1);
            float da_dt = (alpha_f - blades[b].alpha_prev) / dt_si;
            blades[b].alpha_prev = alpha_f;

            blades[b].da_dt_buffer.push_back(da_dt);
            if (blades[b].da_dt_buffer.size() > 10) blades[b].da_dt_buffer.pop_front();
            float da_dt_avg = 0.0f;
            for (float v : blades[b].da_dt_buffer) da_dt_avg += v;
            da_dt_avg /= max(1.0f, (float)blades[b].da_dt_buffer.size());

            float Re_local = units.si_Re(chord_si, units.si_u(V_mag), si_nu);

            // --- DYNAMIC STALL ---
            float alpha_dyn_deg = alpha_f;
            float pitch_term = (chord_si * radians(da_dt_avg)) / (2.0f * units.si_u(V_mag));
            float sqrt_pitch = sqrtf(fabs(pitch_term));
            float ds_shift = ((alpha_f * da_dt_avg > 0) ? 1.0f : 0.0f) * 6.0f * sqrt_pitch;
            alpha_dyn_deg += (da_dt_avg > 0) ? -ds_shift : ds_shift;

            float lookup_a = fabs(alpha_dyn_deg);
            float cl = interpolate_2d(alphas, Re_nums, cls, lookup_a, Re_local);
            float cd = interpolate_2d(alphas, Re_nums, cds, lookup_a, Re_local);
            float cm = interpolate_2d(alphas, Re_nums, cms, lookup_a, Re_local);
            if (alpha_f < 0) { cl = -cl; cm = -cm; }

            // --- ADDED MASS ---
            float da_dt_rad = radians(da_dt_avg);
            float current_alpha_rad = radians(alpha_f);
            float Un_si = units.si_u(V_mag) * sin(current_alpha_rad);
            float dUn_dt = (Un_si - blades[b].Un_prev) / dt_si;
            blades[b].Un_prev = Un_si;

            float V_rel_sq = max(0.0001f, units.si_u(V_mag) * units.si_u(V_mag));
            float C_n_am = -(pi * chord_si * dUn_dt) / (8.0f * V_rel_sq);
            float C_c_am = (pi * chord_si * da_dt_rad * Un_si) / (8.0f * V_rel_sq);

            float sin_a = sin(current_alpha_rad);
            float cos_a = cos(current_alpha_rad);
            cl += C_n_am * cos_a + C_c_am * sin_a;
            cd += C_n_am * sin_a - C_c_am * cos_a;

            // --- 3D Correction ---
            float AR = H_si / chord_si;
            float factor_3d = 1.0f / (1.0f + 2.0f / (AR * 0.9f));
            float cl_3d = cl * factor_3d;
            float cd_3d = cd + (cl_3d * cl_3d) / (pi * AR * 0.9f);

            float dyn_p = 0.5f * V_mag * V_mag * chord_lbm;
            float3 e_d = normalize(rel_vel); 
            float3 e_l = cross(float3(0, 0, 1), e_d); 
            float3 f_blade = dyn_p * (cl_3d * e_l + cd_3d * e_d);

            // Torque usando producto cruz (componente Z)
            float torque_contribution = cross(pos_AC - center, f_blade).z;
            torque_lbm_instant += torque_contribution;

            smear_force(lbm, pos_AC, -f_blade, 0.25f * chord_lbm);

            if (b == 0 && lbm.get_t() >= WARMUP_STEPS && lbm.get_t() % 20 == 0) {
                float Cn = cl_3d * cos_a + cd_3d * sin_a;
                float Ct = cl_3d * sin_a - cd_3d * cos_a;
                blade_file << units.si_t(lbm.get_t()) << ","
                    << degrees(theta) << "," << degrees(phi_norm) << ","
                    << alpha_f << "," << cl_3d << "," << cd_3d << ","
                    << Cn << "," << Ct << "," << units.si_u(V_mag) << ","
                    << units.si_M(torque_contribution) * (H_si / units.si_x(lbm_N.z)) << endl;
            }
        }

        lbm.F.write_to_device();
        lbm.run(1u);

        float theta_old = theta;
        theta = fmod(theta + lbm_omega, 2.0f * pi);
        if (theta_old > pi && theta < pi) torque_samples++; 

        if (lbm.get_t() >= WARMUP_STEPS) {
            float si_torque = units.si_M(torque_lbm_instant) * (H_si / units.si_x(lbm_N.z));
            torque_output_buffer.push_back(si_torque);
            if (torque_output_buffer.size() > 15) torque_output_buffer.pop_front();

            float t_smooth = accumulate(torque_output_buffer.begin(), torque_output_buffer.end(), 0.0f) / torque_output_buffer.size();
            torque_sum += si_torque;
            torque_samples++;

            if (lbm.get_t() % 20 == 0) torque_file << units.si_t(lbm.get_t()) << "\t" << t_smooth << endl;
            if (lbm.get_t() % 100 == 0) {
                double cp = (torque_sum / max(1L, torque_samples)) * si_omega / wind_power_si;
                cout << "Step: " << lbm.get_t() << " | T: " << fixed << setprecision(2) << t_smooth << " | CP: " << setprecision(3) << cp << endl;
            }
        }
        else if (lbm.get_t() % 500 == 0) {
            cout << "WARMUP: " << (100.0f * lbm.get_t() / WARMUP_STEPS) << "%" << endl;
        }
    }

    torque_file.close();
    blade_file.close();
}

[IsTrivial]

I think this is the right approach to avoid the bottleneck I had, but I still have the resolution issue.

#include "setup.hpp"
#include <fstream>
#include <vector>
#include <deque>
#include <numeric>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <CL/opencl.hpp>

using namespace std;

// ============================================================================
// KERNEL OPENCL
// ============================================================================
const string kernel_source = R"(
#define PI 3.14159265359f
#define BUFFER_SIZE 16

typedef struct {
    float alpha_filtered, alpha_prev, Un_prev;
    float da_dt_buffer[BUFFER_SIZE];
    int buffer_idx, buffer_count;
    float cl, cd, cn, ct, v_mag, torque, alpha_deg;
} BladeStateGPU;

int get_index(__global const float* arr, int size, float val) {
    if (val <= arr[0]) return 0;
    if (val >= arr[size-1]) return size - 2;
    for (int i = 0; i < size - 1; i++)
        if (val >= arr[i] && val < arr[i+1]) return i;
    return size - 2;
}

float interpolate_2d_gpu(__global const float* alphas, int n_alpha, __global const float* Res, int n_re,
                         __global const float* table, float alpha, float Re) {
    float a_clim = clamp(alpha, alphas[0], alphas[n_alpha-1]);
    float r_clim = clamp(Re, Res[0], Res[n_re-1]);
    int ai = get_index(alphas, n_alpha, a_clim);
    int ri = get_index(Res, n_re, r_clim);
    float afrac = (a_clim - alphas[ai]) / (alphas[ai+1] - alphas[ai]);
    float Re_frac = (log(r_clim) - log(Res[ri])) / (log(Res[ri+1]) - log(Res[ri]));
    int idx00 = ri * n_alpha + ai, idx01 = idx00 + 1;
    int idx10 = (ri + 1) * n_alpha + ai, idx11 = idx10 + 1;
    float v0 = table[idx00] + afrac * (table[idx01] - table[idx00]);
    float v1 = table[idx10] + afrac * (table[idx11] - table[idx10]);
    return v0 + Re_frac * (v1 - v0);
}

float3 get_vel_soa(__global const float* u_buffer, int grid_size, int Nx, int Ny, float x, float y) {
    int ix = clamp((int)floor(x), 0, Nx - 2), iy = clamp((int)floor(y), 0, Ny - 2);
    float fx = x - ix, fy = y - iy;
    int i00 = iy * Nx + ix, i01 = (iy + 1) * Nx + ix;
    int i10 = iy * Nx + (ix + 1), i11 = (iy + 1) * Nx + (ix + 1);
    int off_y = grid_size, off_z = 2 * grid_size;
    float3 u00 = (float3)(u_buffer[i00], u_buffer[i00+off_y], u_buffer[i00+off_z]);
    float3 u01 = (float3)(u_buffer[i01], u_buffer[i01+off_y], u_buffer[i01+off_z]);
    float3 u10 = (float3)(u_buffer[i10], u_buffer[i10+off_y], u_buffer[i10+off_z]);
    float3 u11 = (float3)(u_buffer[i11], u_buffer[i11+off_y], u_buffer[i11+off_z]);
    return (u00 * (1.0f - fy) + u01 * fy) * (1.0f - fx) + (u10 * (1.0f - fy) + u11 * fy) * fx;
}

__kernel void update_blades(__global float* F_buffer, __global const float* u_buffer, int grid_size,
    __global BladeStateGPU* states, __global const float* alphas, int n_alpha, __global const float* Res, int n_re,
    __global const float* cls, __global const float* cds, __global const float* cms,
    float time_val, float dt_si, float omega, float R_lbm, float chord_lbm, float chord_si,
    float center_x, float center_y, float si_nu, int Nx, int Ny, int Nz, float unit_u, float unit_rho, float unit_x) {
    
    int b = get_global_id(0);
    if (b >= 3) return;
    
    float theta_rotor = fmod(omega * time_val, 2.0f * PI);
    float phi = theta_rotor + 2.0f * PI * b / 3.0f;
    float3 unit_r = (float3)(cos(phi), sin(phi), 0.0f);
    float3 unit_t = (float3)(-sin(phi), cos(phi), 0.0f);
    float3 pos_AC = (float3)(center_x, center_y, 0.5f) + R_lbm * unit_r;
    float3 blade_vel = omega * R_lbm * unit_t;
    
    float3 u_fluid = (float3)(0.0f);
    for(int k=-1; k<=1; k++)
        u_fluid += get_vel_soa(u_buffer, grid_size, Nx, Ny, 
                   (pos_AC + (float)k * 0.25f * chord_lbm * unit_t).x,
                   (pos_AC + (float)k * 0.25f * chord_lbm * unit_t).y);
    u_fluid /= 3.0f;
    
    float3 rel_vel = u_fluid - blade_vel;
    float V_mag = length(rel_vel);
    if (V_mag < 0.001f) return;
    
    float alpha_rad = atan2(dot(rel_vel, unit_r), -dot(rel_vel, unit_t)) + (omega * chord_lbm) / (4.0f * V_mag);
    float alpha_deg = degrees(alpha_rad);
    
    __global BladeStateGPU* s = &states[b];
    if (time_val < 10.0f) s->alpha_filtered = alpha_deg;
    float alpha_f = s->alpha_filtered * 0.92f + alpha_deg * 0.08f;
    s->alpha_filtered = alpha_f;
    float da_dt = (alpha_f - s->alpha_prev) / dt_si;
    s->alpha_prev = alpha_f;
    
    s->da_dt_buffer[s->buffer_idx] = da_dt;
    s->buffer_idx = (s->buffer_idx + 1) % BUFFER_SIZE;
    if (s->buffer_count < BUFFER_SIZE) s->buffer_count++;
    float da_dt_avg = 0.0f;
    for(int i=0; i<s->buffer_count; i++) da_dt_avg += s->da_dt_buffer[i];
    da_dt_avg /= (float)s->buffer_count;
    
    float Re_local = (V_mag * unit_u) * chord_si / si_nu;
    float pitch_term = (chord_si * radians(da_dt_avg)) / (2.0f * (V_mag * unit_u));
    float ds_shift = ((alpha_f * da_dt_avg > 0) ? 1.0f : 0.0f) * 6.0f * sqrt(fabs(pitch_term));
    float alpha_dyn = alpha_f + ((da_dt_avg > 0) ? -ds_shift : ds_shift);
    
    float cl = interpolate_2d_gpu(alphas, n_alpha, Res, n_re, cls, fabs(alpha_dyn), Re_local);
    float cd = interpolate_2d_gpu(alphas, n_alpha, Res, n_re, cds, fabs(alpha_dyn), Re_local);
    if (alpha_f < 0) cl = -cl;
    
    float AR = 1.0f / 0.14f;
    float cl_3d = cl / (1.0f + 2.0f / (AR * 0.9f));
    float cd_3d = cd + (cl_3d * cl_3d) / (PI * AR * 0.9f);
    
    float dyn_p = 0.5f * V_mag * V_mag * chord_lbm;
    float3 e_d = normalize(rel_vel);
    float3 e_l = cross((float3)(0,0,1), e_d);
    float3 f_blade = dyn_p * (cl_3d * e_l + cd_3d * e_d);
    
    s->cl = cl_3d; s->cd = cd_3d;
    s->cn = cl_3d * cos(radians(alpha_f)) + cd_3d * sin(radians(alpha_f));
    s->ct = cl_3d * sin(radians(alpha_f)) - cd_3d * cos(radians(alpha_f));
    s->v_mag = V_mag * unit_u; s->alpha_deg = alpha_f;
    s->torque = cross(pos_AC - (float3)(center_x, center_y, 0.5f), f_blade).z;
    
    float epsilon_lu = 0.25f * chord_lbm, sigma_sq = epsilon_lu * epsilon_lu;
    float gauss_norm = 1.0f / (2.0f * PI * sigma_sq);
    int range = (int)ceil(3.0f * epsilon_lu);
    int min_x = max(0, (int)floor(pos_AC.x) - range), max_x = min(Nx - 1, (int)floor(pos_AC.x) + range);
    int min_y = max(0, (int)floor(pos_AC.y) - range), max_y = min(Ny - 1, (int)floor(pos_AC.y) + range);
    
    for (int xx = min_x; xx <= max_x; ++xx) {
        for (int yy = min_y; yy <= max_y; ++yy) {
            float3 dr = (float3)(xx + 0.5f, yy + 0.5f, 0.5f) - pos_AC;
            float r2 = dot(dr, dr);
            if (r2 > 9.0f * sigma_sq) continue;
            float g = gauss_norm * exp(-r2 / (2.0f * sigma_sq));
            int n = yy * Nx + xx;
            F_buffer[n] -= f_blade.x * g;
            F_buffer[n + grid_size] -= f_blade.y * g;
        }
    }
}
)";

struct BladeStateGPU_Host {
    float alpha_filtered, alpha_prev, Un_prev;
    float da_dt_buffer[16];
    int buffer_idx, buffer_count;
    float cl, cd, cn, ct, v_mag, torque, alpha_deg;
};

// ============================================================================
// CARGA DE DATOS
// ============================================================================
bool read_airfoil_data(const string& filename, vector<float>& alphas, vector<float>& Re_nums,
    vector<vector<float>>& cls, vector<vector<float>>& cds, vector<vector<float>>& cms) {
    ifstream file(filename);
    if (!file) return false;
    string line, token;
    if (!getline(file, line)) return false;
    stringstream ss(line);
    getline(ss, token, ',');
    while (getline(ss, token, ',')) {
        token.erase(token.find_last_not_of(" \n\r\t") + 1);
        if (token == "|" || token.empty()) break;
        alphas.push_back(stof(token));
    }
    while (getline(file, line)) {
        stringstream row(line);
        getline(row, token, ',');
        if (token.empty()) continue;
        try { Re_nums.push_back(stof(token)); }
        catch (...) { continue; }
        vector<float> cl_row, cd_row, cm_row;
        int section = 0;
        while (getline(row, token, ',')) {
            token.erase(token.find_last_not_of(" \n\r\t") + 1);
            if (token == "|") { section++; continue; }
            if (token.empty()) continue;
            float val = stof(token);
            if (section == 0) cl_row.push_back(val);
            else if (section == 1) cd_row.push_back(val);
            else if (section == 2) cm_row.push_back(val);
        }
        cls.push_back(cl_row); cds.push_back(cd_row);
        cms.push_back(cm_row.empty() ? vector<float>(cl_row.size(), 0.0f) : cm_row);
    }
    return !Re_nums.empty();
}

// ============================================================================
// MAIN SETUP
// ============================================================================
void main_setup() {
    cout << unitbuf; cerr << unitbuf;
    cout << "=== INICIALIZANDO PARAMETROS (GPU ACELERADO) ===" << endl;

    const float D_si = 1.0f, H_si = 1.0f, chord_si = 0.14f, U_inf = 1.0f, rpm = 36.3f;
    const int N_blades = 3;
    const float si_nu = 1.0e-6f, si_rho = 998.0f, R_si = D_si / 2.0f;
    const float si_omega = 2.0f * pi * (rpm / 60.0f), si_period = 2.0f * pi / si_omega;
    const float domain_x_si = 35.0f * D_si, domain_y_si = 20.0f * D_si;
    const uint3 lbm_N = resolution(float3(domain_x_si, domain_y_si, 1.0f), 1500u);
    const float lbm_u = 0.1f;

    units.set_m_kg_s((float)lbm_N.x, lbm_u, 1.0f, domain_x_si, U_inf, si_rho);
    const float lbm_nu = units.nu(si_nu), lbm_R = units.x(R_si), chord_lbm = units.x(chord_si);
    const ulong lbm_T = units.t(15.0f * si_period), WARMUP_STEPS = 7000;
    const float lbm_omega = units.omega(si_omega);

    LBM lbm(lbm_N, 1u, 1u, 1u, lbm_nu);

    parallel_for(lbm.get_N(), [&](ulong n) {
        lbm.rho[n] = 1.0f; lbm.u.x[n] = lbm_u; lbm.u.y[n] = 0.0f; lbm.u.z[n] = 0.0f;
        uint x, y, z; lbm.coordinates(n, x, y, z);
        if (x == 0u) {
            lbm.flags[n] = TYPE_E;
            lbm.u.x[n] = lbm_u; 
        }
        else if (x == lbm.get_Nx() - 1u) {
            lbm.flags[n] = TYPE_E;
        }
        });

    lbm.rho.write_to_device(); lbm.u.write_to_device(); lbm.flags.write_to_device();

    vector<float> alphas, Re_nums;
    vector<vector<float>> cls, cds, cms;
    if (!read_airfoil_data("airfoil_data.csv", alphas, Re_nums, cls, cds, cms)) return;

    vector<float> cls_flat, cds_flat, cms_flat;
    for (const auto& row : cls) cls_flat.insert(cls_flat.end(), row.begin(), row.end());
    for (const auto& row : cds) cds_flat.insert(cds_flat.end(), row.begin(), row.end());
    for (const auto& row : cms) cms_flat.insert(cms_flat.end(), row.begin(), row.end());

    // OpenCL Setup
    const cl::Buffer& buf_u_fluid = lbm.lbm_domain[0]->u.get_cl_buffer();
    cl::Context context = buf_u_fluid.getInfo<CL_MEM_CONTEXT>();
    cl::Device device = context.getInfo<CL_CONTEXT_DEVICES>()[0];
    cl::CommandQueue queue(context, device);

    cl::Program::Sources sources; sources.push_back({ kernel_source.c_str(), kernel_source.length() });
    cl::Program program(context, sources);
    if (program.build({ device }) != CL_SUCCESS) {
        cout << "Error building: " << program.getBuildInfo<CL_PROGRAM_BUILD_LOG>(device) << endl;
        return;
    }
    cl::Kernel blade_kernel_raw(program, "update_blades");

    cl::Buffer buf_alphas(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * alphas.size(), alphas.data());
    cl::Buffer buf_Re(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * Re_nums.size(), Re_nums.data());
    cl::Buffer buf_cls(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * cls_flat.size(), cls_flat.data());
    cl::Buffer buf_cds(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * cds_flat.size(), cds_flat.data());
    cl::Buffer buf_cms(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * cms_flat.size(), cms_flat.data());

    vector<BladeStateGPU_Host> host_states(N_blades, { 0 });
    cl::Buffer buf_states(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, sizeof(BladeStateGPU_Host) * N_blades, host_states.data());

    const cl::Buffer& buf_F_fluid_const = lbm.lbm_domain[0]->F.get_cl_buffer();
    cl::Buffer& buf_F_fluid = const_cast<cl::Buffer&>(buf_F_fluid_const);

    int grid_size = (int)(lbm.get_Nx() * lbm.get_Ny() * lbm.get_Nz());
    const float center_x = units.x(10.0f * D_si), center_y = lbm.center().y;

    ofstream torque_file("torque_data.txt"), blade_file("blade_debug.csv");
    torque_file << "# Time\tTorque_Smooth" << endl;
    blade_file << "Time_s,Theta_deg,Phi_deg,Alpha_deg,Cl_3d,Cd_3d,Cn,Ct,V_mag,Torque_inst" << endl;

    deque<float> torque_output_buffer;
    double torque_sum = 0.0; long torque_samples = 0;
    const float wind_power_si = 0.5f * si_rho * (D_si * H_si) * pow(U_inf, 3.0f);

    cout << "=== STARTING MAIN LOOP ===" << endl;
    lbm.run(0u);

    while (lbm.get_t() < lbm_T) {
        float zero = 0.0f;
        queue.enqueueFillBuffer(buf_F_fluid, zero, 0, sizeof(float) * grid_size * 3);

        blade_kernel_raw.setArg(0, buf_F_fluid); blade_kernel_raw.setArg(1, buf_u_fluid);
        blade_kernel_raw.setArg(2, grid_size); blade_kernel_raw.setArg(3, buf_states);
        blade_kernel_raw.setArg(4, buf_alphas); blade_kernel_raw.setArg(5, (int)alphas.size());
        blade_kernel_raw.setArg(6, buf_Re); blade_kernel_raw.setArg(7, (int)Re_nums.size());
        blade_kernel_raw.setArg(8, buf_cls); blade_kernel_raw.setArg(9, buf_cds);
        blade_kernel_raw.setArg(10, buf_cms); blade_kernel_raw.setArg(11, (float)lbm.get_t());
        blade_kernel_raw.setArg(12, units.si_t(1)); blade_kernel_raw.setArg(13, lbm_omega);
        blade_kernel_raw.setArg(14, lbm_R); blade_kernel_raw.setArg(15, chord_lbm);
        blade_kernel_raw.setArg(16, chord_si); blade_kernel_raw.setArg(17, center_x);
        blade_kernel_raw.setArg(18, center_y); blade_kernel_raw.setArg(19, si_nu);
        blade_kernel_raw.setArg(20, (int)lbm.get_Nx()); blade_kernel_raw.setArg(21, (int)lbm.get_Ny());
        blade_kernel_raw.setArg(22, (int)lbm.get_Nz()); blade_kernel_raw.setArg(23, units.si_u(1.0f));
        blade_kernel_raw.setArg(24, units.si_rho(1.0f)); blade_kernel_raw.setArg(25, units.si_x(1.0f));

        queue.enqueueNDRangeKernel(blade_kernel_raw, cl::NullRange, cl::NDRange(N_blades), cl::NullRange);
        queue.finish();
        lbm.run(1u);

        if (lbm.get_t() >= WARMUP_STEPS || lbm.get_t() % 20 == 0) {
            queue.enqueueReadBuffer(buf_states, CL_TRUE, 0, sizeof(BladeStateGPU_Host) * N_blades, host_states.data());

            float total_torque_lbm = 0.0f;
            for (int b = 0; b < N_blades; b++) total_torque_lbm += host_states[b].torque;
            float si_torque = units.si_M(total_torque_lbm) * (H_si / units.si_x(lbm_N.z));

            if (lbm.get_t() >= WARMUP_STEPS) {
                torque_output_buffer.push_back(si_torque);
                if (torque_output_buffer.size() > 15) torque_output_buffer.pop_front();
                float t_smooth = accumulate(torque_output_buffer.begin(), torque_output_buffer.end(), 0.0f) / torque_output_buffer.size();
                torque_sum += si_torque; torque_samples++;

                if (lbm.get_t() % 20 == 0) torque_file << units.si_t(lbm.get_t()) << "\t" << t_smooth << endl;
                if (lbm.get_t() % 100 == 0) {
                    double cp = (torque_sum / max(1L, torque_samples)) * si_omega / wind_power_si;
                    cout << "Step: " << lbm.get_t() << " | T: " << fixed << setprecision(2) << t_smooth
                        << " | CP: " << setprecision(3) << cp << endl;
                }
            }
            else if (lbm.get_t() % 500 == 0) {
                cout << "WARMUP: " << (100.0f * lbm.get_t() / WARMUP_STEPS) << "%" << endl;
            }

            if (lbm.get_t() % 20 == 0) {
                float theta_rad = fmod(lbm_omega * (float)lbm.get_t(), 2.0f * pi);
                float theta_deg = theta_rad * (180.0f / pi);
                blade_file << units.si_t(lbm.get_t()) << "," << theta_deg << "," << theta_deg << ","
                    << host_states[0].alpha_deg << "," << host_states[0].cl << "," << host_states[0].cd << ","
                    << host_states[0].cn << "," << host_states[0].ct << "," << host_states[0].v_mag << ","
                    << units.si_M(host_states[0].torque) * (H_si / units.si_x(lbm_N.z)) << endl;
            }
        }
    }

    torque_file.close(); blade_file.close();
}

[jasonxauat]

I think ALM and IBM are quite similar in implementation. You can refer to the IBM code I uploaded, which can already calculate simple rigid body motion in GPU kernel code. However, the calculation of total forces and moments remains performed on the CPU, with only a small amount of data needing to be transferred at each time step. #306