Help with the implementation of this ALM over a H-VAWT Darrieus in water

[IsTrivial]

I’ve been implementing an Actuator Line Model (ALM) for an H-VAWT in water, following a published paper, to evaluate how accurately this method predicts output torque compared to experimental results. To build the model, I collected data for multiple angles of attack and Reynolds numbers and interpolated the polar coefficients.

My current results are:

Time 49700: ALM Torque = 16.7655 N·m, Avg ALM Torque = 44.9571 N·m, Avg Power = 282.4738 W, Wake Torque = 41.9914 N·m
Wake Velocity at 3 radius (~1.35m downstream): 1.4412 m/s (expected <1.62 m/s)
Wake Velocity at 5 radius (~2.25m downstream): 1.4337 m/s (expected deficit, recovering)
Wake Velocity at 10 radius (~4.5m downstream): 1.4535 m/s (expected near 1.62 m/s)
Average domain velocity (x-component, upstream only): 1.6958 m/s (expected ~1.62 m/s in steady state)

I’m trying to determine whether there are major flaws in my implementation, or if I might be misunderstanding certain FluidX3D functions used in my code. Any feedback or corrections would be greatly appreciated.

#include "defines.hpp"
#include "lbm.hpp"
#include "shapes.hpp"
#include <cmath>
#include <cstddef>
#include <algorithm>
#include <vector>
#include <fstream>
#include <chrono>
#include <deque>
#include <numeric>
#include <tuple>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// 1D interpolation function
float interpolate(const std::vector<float>& x, const std::vector<float>& y, float val) {
    if (val <= x.front()) return y.front();
    if (val >= x.back()) return y.back();
    auto it = std::lower_bound(x.begin(), x.end(), val);
    size_t i = it - x.begin() - 1;
    float frac = (val - x[i]) / (x[i + 1] - x[i]);
    return y[i] + frac * (y[i + 1] - y[i]);
}

// 2D interpolation function for cl and cd
float interpolate_2d(const std::vector<float>& alphas, const std::vector<float>& Re_nums, 
                    const std::vector<std::vector<float>>& values, float alpha, float Re) {
    // Clamp alpha and Re to table bounds
    float alpha_clamped = std::max(alphas.front(), std::min(alphas.back(), alpha));
    float Re_clamped = std::max(Re_nums.front(), std::min(Re_nums.back(), Re));

    // Find indices for alpha
    auto alpha_it = std::lower_bound(alphas.begin(), alphas.end(), alpha_clamped);
    size_t alpha_i = alpha_it - alphas.begin() - 1;
    if (alpha_i >= alphas.size() - 1) alpha_i = alphas.size() - 2;
    float alpha_frac = (alpha_clamped - alphas[alpha_i]) / (alphas[alpha_i + 1] - alphas[alpha_i]);

    // Find indices for Re
    auto Re_it = std::lower_bound(Re_nums.begin(), Re_nums.end(), Re_clamped);
    size_t Re_i = Re_it - Re_nums.begin() - 1;
    if (Re_i >= Re_nums.size() - 1) Re_i = Re_nums.size() - 2;
    float Re_frac = (Re_clamped - Re_nums[Re_i]) / (Re_nums[Re_i + 1] - Re_nums[Re_i]);

    // Bilinear interpolation
    float v00 = values[Re_i][alpha_i];
    float v01 = values[Re_i][alpha_i + 1];
    float v10 = values[Re_i + 1][alpha_i];
    float v11 = values[Re_i + 1][alpha_i + 1];

    float v0 = v00 + alpha_frac * (v01 - v00);
    float v1 = v10 + alpha_frac * (v11 - v10);
    return v0 + Re_frac * (v1 - v0);
}

float3 get_interpolated_velocity(LBM& lbm, float3 pos) {
    int ix = (int)std::floor(pos.x);
    int iy = (int)std::floor(pos.y);
    int iz = (int)std::floor(pos.z);
    float fx = pos.x - ix;
    float fy = pos.y - iy;
    float fz = pos.z - iz;

    // Clamp to boundaries
    ix = std::max(0, std::min(ix, (int)lbm.get_Nx() - 2));
    iy = std::max(0, std::min(iy, (int)lbm.get_Ny() - 2));
    iz = std::max(0, std::min(iz, (int)lbm.get_Nz() - 2));

    float3 u000 = float3(lbm.u[lbm.index(ix, iy, iz)], lbm.u[lbm.index(ix, iy, iz) + lbm.get_N()], lbm.u[lbm.index(ix, iy, iz) + 2 * lbm.get_N()]);
    float3 u001 = float3(lbm.u[lbm.index(ix, iy, iz + 1)], lbm.u[lbm.index(ix, iy, iz + 1) + lbm.get_N()], lbm.u[lbm.index(ix, iy, iz + 1) + 2 * lbm.get_N()]);
    float3 u010 = float3(lbm.u[lbm.index(ix, iy + 1, iz)], lbm.u[lbm.index(ix, iy + 1, iz) + lbm.get_N()], lbm.u[lbm.index(ix, iy + 1, iz) + 2 * lbm.get_N()]);
    float3 u011 = float3(lbm.u[lbm.index(ix, iy + 1, iz + 1)], lbm.u[lbm.index(ix, iy + 1, iz + 1) + lbm.get_N()], lbm.u[lbm.index(ix, iy + 1, iz + 1) + 2 * lbm.get_N()]);
    float3 u100 = float3(lbm.u[lbm.index(ix + 1, iy, iz)], lbm.u[lbm.index(ix + 1, iy, iz) + lbm.get_N()], lbm.u[lbm.index(ix + 1, iy, iz) + 2 * lbm.get_N()]);
    float3 u101 = float3(lbm.u[lbm.index(ix + 1, iy, iz + 1)], lbm.u[lbm.index(ix + 1, iy, iz + 1) + lbm.get_N()], lbm.u[lbm.index(ix + 1, iy, iz + 1) + 2 * lbm.get_N()]);
    float3 u110 = float3(lbm.u[lbm.index(ix + 1, iy + 1, iz)], lbm.u[lbm.index(ix + 1, iy + 1, iz) + lbm.get_N()], lbm.u[lbm.index(ix + 1, iy + 1, iz) + 2 * lbm.get_N()]);
    float3 u111 = float3(lbm.u[lbm.index(ix + 1, iy + 1, iz + 1)], lbm.u[lbm.index(ix + 1, iy + 1, iz + 1) + lbm.get_N()], lbm.u[lbm.index(ix + 1, iy + 1, iz + 1) + 2 * lbm.get_N()]);

    float3 u00 = u000 * (1 - fz) + u001 * fz;
    float3 u01 = u010 * (1 - fz) + u011 * fz;
    float3 u10 = u100 * (1 - fz) + u101 * fz;
    float3 u11 = u110 * (1 - fz) + u111 * fz;
    float3 u0 = u00 * (1 - fy) + u01 * fy;
    float3 u1 = u10 * (1 - fy) + u11 * fy;
    return u0 * (1 - fx) + u1 * fx;
}

struct MLDataPoint {
    float time;
    float theta;
    float torque;
    float TSR;
    float Re_avg;
    float inflow_velocity;
    float wake_deficit;
};

void main_setup() {
    // ==================== SIMULATION SETUP ====================
    const uint3 lbm_N = resolution(float3(80.0f, 32.0f, 8.0f), 3000u);
    const float si_length = 14.4f * 2.0f;
    const float si_velocity = 1.62f;
    const float si_density = 1000.0f;
    const float si_nu = 1.0e-6f;
    const float lbm_u = 0.04f;
    units.set_m_kg_s((float)lbm_N.x, lbm_u, 1.0f, si_length, si_velocity, si_density);
    const float lbm_nu = units.nu(si_nu);
    const float dt = units.si_t(1.0f);
    const float si_rpm = 60.0f;
    const float si_omega = 2.0f * M_PI * (si_rpm / 60.0f);
    const float radius_si = 0.45f;
    const float lbm_radius = radius_si / units.si_x(1.0f);
    const float lbm_omega = -si_omega * dt;
    const ulong lbm_T = 120000ull;
    const ulong lbm_dt = 5ull;

    printf("=== Conversion Factors ===\n");
    printf("SI length per cell: %.6f m\n", units.si_x(1.0f));
    printf("SI time per step: %.6f s\n", dt);
    printf("SI velocity per LBM u: %.6f m/s\n", units.si_u(1.0f));
    printf("LBM velocity should be: %.6f (SI: %.6f m/s)\n", lbm_u, units.si_u(lbm_u));
    printf("LBM viscosity: %.6e, SI viscosity: %.6e m²/s\n", lbm_nu, units.si_nu(lbm_nu));
    printf("Turbine diameter in LBM units: %.4f, SI units: %.4f m\n", 2.0f * lbm_radius, units.si_x(2.0f * lbm_radius));

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

    // ==================== GEOMETRY SETUP ====================
    const float3 center = float3(12.0f * lbm_radius, lbm.center().y, lbm.center().z);
    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) {
            lbm.flags[n] = TYPE_E;
            lbm.u[n] = lbm_u;
            lbm.u[n + lbm.get_N()] = 0.0f;
            lbm.u[n + 2 * lbm.get_N()] = 0.0f;
        } else if (x == Nx - 1u) {
            lbm.flags[n] = TYPE_E;
            lbm.rho[n] = 1.0f;
        } else if (y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) {
            lbm.flags[n] = TYPE_E;
            lbm.rho[n] = 1.0f;
            lbm.u[n] = lbm_u;
            lbm.u[n + lbm.get_N()] = 0.0f;
            lbm.u[n + 2 * lbm.get_N()] = 0.0f;
        }
    });

    parallel_for(lbm.get_N(), [&](ulong n) {
        lbm.rho[n] = 1.0f;
        lbm.u[n] = lbm_u;
        lbm.u[n + lbm.get_N()] = 0.0f;
        lbm.u[n + 2 * lbm.get_N()] = 0.0f;
    });

    // ==================== DARRIEUS VAWT ACTUATOR LINE SETUP ====================
    const int num_blades = 3;
    const float chord_si = 0.13275f;
    const float chord_lbm = chord_si / units.si_x(1.0f);
    const float blade_height_si = 0.7f;
    const float lbm_blade_height = blade_height_si / units.si_x(1.0f);
    const float z_start = center.z - 0.5f * lbm_blade_height;
    const int num_segments = 100;
    const float dz = lbm_blade_height / (float)num_segments;
    const float sigma = 0.25f * chord_lbm;

    // NACA 0025 airfoil data
    std::vector<float> alphas = {
        0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f, 17.0f, 18.0f, 19.0f, 20.0f, 21.0f, 22.0f, 23.0f, 24.0f, 25.0f, 26.0f, 27.0f, 30.0f, 35.0f, 40.0f, 45.0f, 50.0f, 55.0f, 60.0f, 65.0f, 70.0f, 75.0f, 80.0f, 85.0f, 90.0f, 95.0f, 100.0f, 105.0f, 110.0f, 115.0f, 120.0f, 125.0f, 130.0f, 135.0f, 140.0f, 145.0f, 150.0f, 155.0f, 160.0f, 165.0f, 170.0f, 175.0f, 180.0f
    };
    std::vector<float> Re_nums = {
        10000.0f, 20000.0f, 40000.0f, 80000.0f, 160000.0f, 360000.0f, 700000.0f, 1000000.0f, 2000000.0f, 5000000.0f
    };
    std::vector<std::vector<float>> cls = {
        {0.0000f, -0.0434f, -0.0807f, -0.1151f, -0.1403f, -0.1585f, -0.1723f, -0.1799f, -0.1797f, -0.1776f, -0.1747f, -0.1647f, -0.1512f, -0.1230f, -0.1134f, -0.0604f, 0.0082f, 0.0883f, 0.1785f, 0.3205f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f},
        {0.0000f, -0.0223f, -0.0376f, -0.0509f, -0.0550f, -0.0589f, -0.0650f, -0.0690f, -0.0683f, -0.0653f, -0.0635f, -0.0607f, -0.0572f, -0.0377f, -0.0027f, 0.0512f, 0.1213f, 0.2071f, 0.3435f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f},
        {0.0000f, 0.0411f, 0.0821f, 0.1167f, 0.1479f, 0.1693f, 0.1778f, 0.1856f, 0.1929f, 0.1942f, 0.1915f, 0.1889f, 0.1813f, 0.1753f, 0.1669f, 0.1824f, 0.2183f, 0.2843f, 0.4026f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f},
        {0.0000f, 0.0743f, 0.1487f, 0.2198f, 0.2832f, 0.3413f, 0.3866f, 0.4215f, 0.4503f, 0.4656f, 0.4749f, 0.4797f, 0.4789f, 0.4754f, 0.4713f, 0.4588f, 0.4522f, 0.4547f, 0.4826f, 0.5574f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f},
        {0.0000f, 0.0875f, 0.1745f, 0.2613f, 0.3400f, 0.4170f, 0.4831f, 0.5263f, 0.5916f, 0.6271f, 0.6539f, 0.6720f, 0.6811f, 0.6866f, 0.6886f, 0.6906f, 0.6921f, 0.6935f, 0.6945f, 0.7006f, 0.7069f, 0.7164f, 0.7289f, 0.7421f, 0.7605f, 0.7840f, 0.8047f, 0.8306f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f},
        {0.0000f, 0.0982f, 0.1721f, 0.2682f, 0.3608f, 0.4495f, 0.5295f, 0.6062f, 0.6701f, 0.7256f, 0.7710f, 0.8039f, 0.8295f, 0.8471f, 0.8565f, 0.8656f, 0.8742f, 0.8821f, 0.8890f, 0.8992f, 0.9114f, 0.9216f, 0.9381f, 0.9529f, 0.9721f, 0.9907f, 1.0147f, 1.0396f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f},
        {0.0000f, 0.0817f, 0.1851f, 0.2861f, 0.3811f, 0.4743f, 0.5627f, 0.6440f, 0.7194f, 0.7852f, 0.8405f, 0.8873f, 0.9230f, 0.9491f, 0.9685f, 0.9820f, 0.9944f, 1.0052f, 1.0157f, 1.0276f, 1.0403f, 1.0521f, 1.0678f, 1.0831f, 1.1035f, 1.1223f, 1.1458f, 1.1704f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f},
        {0.0000f, 0.1100f, 0.1890f, 0.2917f, 0.3900f, 0.4844f, 0.5757f, 0.6622f, 0.7387f, 0.8117f, 0.8710f, 0.9228f, 0.9656f, 0.9960f, 1.0192f, 1.0382f, 1.0522f, 1.0654f, 1.0774f, 1.0896f, 1.1031f, 1.1151f, 1.1294f, 1.1476f, 1.1645f, 1.1856f, 1.2071f, 1.2317f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f},
        {0.0000f, 0.1100f, 0.2200f, 0.2978f, 0.4000f, 0.4911f, 0.5913f, 0.6810f, 0.7664f, 0.8415f, 0.9135f, 0.9723f, 1.0245f, 1.0672f, 1.0989f, 1.1258f, 1.1488f, 1.1654f, 1.1806f, 1.1953f, 1.2094f, 1.2220f, 1.2354f, 1.2537f, 1.2683f, 1.2890f, 1.3087f, 1.3347f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f},
        {0.0000f, 0.1100f, 0.2200f, 0.3300f, 0.4078f, 0.5097f, 0.6076f, 0.6997f, 0.7886f, 0.8730f, 0.9472f, 1.0177f, 1.0750f, 1.1270f, 1.1717f, 1.2065f, 1.2358f, 1.2617f, 1.2790f, 1.2952f, 1.3121f, 1.3257f, 1.3378f, 1.3558f, 1.3735f, 1.3960f, 1.4187f, 1.4395f, 0.8550f, 0.9800f, 1.0350f, 1.0500f, 1.0200f, 0.9550f, 0.8750f, 0.7600f, 0.6300f, 0.5000f, 0.3650f, 0.2300f, 0.0900f, -0.0500f, -0.1850f, -0.3200f, -0.4500f, -0.5750f, -0.6700f, -0.7600f, -0.8500f, -0.9300f, -0.9800f, -0.9000f, -0.7700f, -0.6700f, -0.6350f, -0.6800f, -0.8500f, -0.6600f, 0.0000f}
    };
    std::vector<std::vector<float>> cds = {
        {0.0455f, 0.0457f, 0.0463f, 0.0474f, 0.0488f, 0.0508f, 0.0532f, 0.0562f, 0.0597f, 0.0638f, 0.0686f, 0.0820f, 0.1230f, 0.1580f, 0.1960f, 0.2380f, 0.2820f, 0.3290f, 0.4050f, 0.5700f, 0.7450f, 0.9200f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f},
        {0.0341f, 0.0342f, 0.0346f, 0.0354f, 0.0365f, 0.0379f, 0.0397f, 0.0419f, 0.0445f, 0.0476f, 0.0511f, 0.0810f, 0.1230f, 0.1580f, 0.1960f, 0.2380f, 0.2820f, 0.3290f, 0.4050f, 0.5700f, 0.7450f, 0.9200f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f},
        {0.0258f, 0.0259f, 0.0263f, 0.0269f, 0.0278f, 0.0290f, 0.0304f, 0.0322f, 0.0343f, 0.0367f, 0.0396f, 0.0428f, 0.0850f, 0.1580f, 0.1960f, 0.2380f, 0.2820f, 0.3290f, 0.4050f, 0.5700f, 0.7450f, 0.9200f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f},
        {0.0197f, 0.0198f, 0.0201f, 0.0207f, 0.0214f, 0.0224f, 0.0236f, 0.0250f, 0.0268f, 0.0288f, 0.0311f, 0.0338f, 0.0368f, 0.0401f, 0.1010f, 0.1960f, 0.2380f, 0.2820f, 0.3290f, 0.4050f, 0.5700f, 0.7450f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f},
        {0.0155f, 0.0156f, 0.0158f, 0.0163f, 0.0169f, 0.0177f, 0.0187f, 0.0201f, 0.0213f, 0.0230f, 0.0250f, 0.0272f, 0.0296f, 0.0324f, 0.0355f, 0.0389f, 0.0426f, 0.0466f, 0.0510f, 0.0558f, 0.0610f, 0.0668f, 0.0730f, 0.0796f, 0.0870f, 0.0947f, 0.1029f, 0.1117f, 0.5700f, 0.7450f, 0.9200f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f},
        {0.0121f, 0.0122f, 0.0125f, 0.0129f, 0.0135f, 0.0141f, 0.0149f, 0.0159f, 0.0170f, 0.0185f, 0.0200f, 0.0219f, 0.0239f, 0.0261f, 0.0286f, 0.0314f, 0.0343f, 0.0375f, 0.0411f, 0.0448f, 0.0489f, 0.0534f, 0.0582f, 0.0634f, 0.0690f, 0.0750f, 0.0814f, 0.0881f, 0.5700f, 0.7450f, 0.9200f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f},
        {0.0104f, 0.0104f, 0.0105f, 0.0109f, 0.0113f, 0.0119f, 0.0126f, 0.0134f, 0.0144f, 0.0156f, 0.0170f, 0.0185f, 0.0202f, 0.0222f, 0.0244f, 0.0267f, 0.0292f, 0.0319f, 0.0349f, 0.0381f, 0.0415f, 0.0453f, 0.0493f, 0.0536f, 0.0583f, 0.0632f, 0.0685f, 0.0741f, 0.5700f, 0.7450f, 0.9200f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f},
        {0.0097f, 0.0098f, 0.0098f, 0.0101f, 0.0105f, 0.0110f, 0.0116f, 0.0124f, 0.0133f, 0.0144f, 0.0157f, 0.0171f, 0.0187f, 0.0205f, 0.0225f, 0.0246f, 0.0269f, 0.0294f, 0.0322f, 0.0351f, 0.0382f, 0.0417f, 0.0454f, 0.0493f, 0.0535f, 0.0581f, 0.0629f, 0.0680f, 0.5700f, 0.7450f, 0.9200f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f},
        {0.0090f, 0.0090f, 0.0091f, 0.0093f, 0.0096f, 0.0100f, 0.0106f, 0.0112f, 0.0120f, 0.0129f, 0.0140f, 0.0151f, 0.0165f, 0.0180f, 0.0197f, 0.0216f, 0.0235f, 0.0257f, 0.0280f, 0.0305f, 0.0332f, 0.0362f, 0.0394f, 0.0428f, 0.0465f, 0.0504f, 0.0545f, 0.0589f, 0.5700f, 0.7450f, 0.9200f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f},
        {0.0084f, 0.0084f, 0.0085f, 0.0087f, 0.0089f, 0.0092f, 0.0097f, 0.0103f, 0.0110f, 0.0117f, 0.0126f, 0.0137f, 0.0148f, 0.0161f, 0.0175f, 0.0190f, 0.0207f, 0.0225f, 0.0245f, 0.0267f, 0.0290f, 0.0316f, 0.0344f, 0.0374f, 0.0404f, 0.0436f, 0.0470f, 0.0509f, 0.5700f, 0.7450f, 0.9200f, 1.0750f, 1.2150f, 1.3450f, 1.4700f, 1.5750f, 1.6650f, 1.7350f, 1.7800f, 1.8000f, 1.8000f, 1.7800f, 1.7500f, 1.7000f, 1.6350f, 1.5550f, 1.4650f, 1.3500f, 1.2250f, 1.0850f, 0.9250f, 0.7550f, 0.5750f, 0.4200f, 0.3200f, 0.2300f, 0.1400f, 0.0550f, 0.0250f}
    };

    // ==================== MAIN SIMULATION LOOP ====================
    lbm.graphics.visualization_modes = VIS_FLAG_LATTICE | VIS_Q_CRITERION;
    lbm.run(0u);
    std::ofstream torque_file("torque_data.txt");
    torque_file << "# Time\tTorque\n";

    const uint x_measure_3 = (uint)(center.x + 3.0f * lbm_radius);
    const uint x_measure_5 = (uint)(center.x + 5.0f * lbm_radius);
    const uint x_measure_10 = (uint)(center.x + 10.0f * lbm_radius);
    float torque_sum = 0.0f;
    ulong measurement_count = 0ull;
    const float dx = units.si_x(1.0f);
    const float cell_area = dx * dx;
    const float swept_area = 0.9f * 0.7f;
    float theta = 0.0f;
    std::deque<float> rev_avg_torque_history;
    const size_t history_size = 10;
    const float convergence_threshold = 0.005f;
    bool converged = false;
    std::vector<float> torque_in_rev;
    float last_theta = 0.0f;
    const float rev_angle = 2.0f * M_PI;
    std::vector<MLDataPoint> ml_dataset;
    auto start_time = std::chrono::high_resolution_clock::now();

    while (lbm.get_t() < lbm_T && !converged) {
    // Read velocity field from device
    lbm.u.read_from_device();
    
    parallel_for(lbm.get_N(), [&](ulong n) {
        uint x = 0u, y = 0u, z = 0u; 
        lbm.coordinates(n, x, y, z);
        
        if (x == Nx - 1u) { // Outlet boundary
            ulong n_interior = lbm.index(Nx - 2u, y, z);
            lbm.u[n] = lbm.u[n_interior];
            lbm.u[n + lbm.get_N()] = lbm.u[n_interior + lbm.get_N()];
            lbm.u[n + 2 * lbm.get_N()] = lbm.u[n_interior + 2 * lbm.get_N()];
            lbm.rho[n] = 1.0f;
        }
    });
    
    lbm.u.write_to_device();

    // Clear force field
    parallel_for(lbm.get_N(), [&](ulong n) {
        lbm.F[n] = 0.0f;
        lbm.F[n + lbm.get_N()] = 0.0f;
        lbm.F[n + 2 * lbm.get_N()] = 0.0f;
    });
    


        float torque_lbm = 0.0f;
        for (int b = 0; b < num_blades; ++b) {
            float phi = theta + 2.0f * M_PI * (float)b / (float)num_blades;
            float cos_phi = cos(phi);
            float sin_phi = sin(phi);
            float3 unit_radial = float3(cos_phi, sin_phi, 0.0f);
            float3 unit_tan = float3(-sin_phi, cos_phi, 0.0f);
            float3 unit_normal = -unit_radial;

            for (int s = 0; s < num_segments; ++s) {
                float z_pos = z_start + (float)s * dz + 0.5f * dz;
                float3 pos = center + float3(lbm_radius * cos_phi, lbm_radius * sin_phi, z_pos - center.z);
                float3 blade_vel = cross(float3(0.0f, 0.0f, lbm_omega), pos - center);
                float3 u_lbm = get_interpolated_velocity(lbm, pos);
                float3 rel_vel = u_lbm - blade_vel;
                float V_t = dot(rel_vel, unit_tan);
                float V_n = dot(rel_vel, unit_normal);
                float V_mag = sqrt(V_t * V_t + V_n * V_n);
                if (V_mag < 1e-6f) continue;

                float alpha_rad;
                if (fabs(V_t) < 1e-6f) {
                    alpha_rad = (V_n > 0.0f) ? M_PI / 2.0f : -M_PI / 2.0f;
                } else {
                    alpha_rad = atan2(V_n, V_t);
                }
                float alpha_deg = alpha_rad * 180.0f / M_PI;
                float abs_alpha = fabs(alpha_deg);

                // Interpolate cl and cd based on alpha and Re
                float Re_local_si = units.si_u(V_mag) * chord_si / si_nu;
                float cl = interpolate_2d(alphas, Re_nums, cls, abs_alpha, Re_local_si);
                float cd = interpolate_2d(alphas, Re_nums, cds, abs_alpha, Re_local_si);
                if (alpha_deg < 0.0f) cl = -cl;

                // VAWT vertical tip loss correction
                float span_fraction = fabsf(z_pos - center.z) / (0.5f * lbm_blade_height);
                float F_tip_vertical = 1.0f;
                if (span_fraction > 0.7f) {
                    float B = (float)num_blades;
                    float arg = exp(-B * (1.0f - span_fraction) / (2.0f * span_fraction));
                    float F_prandtl = (2.0f / M_PI) * static_cast<float>(acos(arg));
                    F_tip_vertical = std::min(1.0f, std::max(0.3f, F_prandtl));
                }
                cl *= F_tip_vertical;
                cd *= F_tip_vertical;

                // Force coefficients
                float cos_alpha = static_cast<float>(cos(alpha_rad));
                float sin_alpha = static_cast<float>(sin(alpha_rad));
                float cn = cl * cos_alpha + cd * sin_alpha;
                float ct = cl * sin_alpha - cd * cos_alpha;

                float dyn_p_times_area = 0.5f * 1.0f * V_mag * V_mag * chord_lbm * dz;

                // VAWT 3D correction
                float aspect_ratio = lbm_blade_height / (2.0f * lbm_radius);
                float TSR_local = (lbm_radius * fabsf(lbm_omega)) / V_mag;
                float solidity = (num_blades * chord_lbm) / (2.0f * M_PI * lbm_radius);
                float temp = cl * cos_alpha + cd * sin_alpha;
                float induction_factor = std::min(0.4f, solidity * temp / 4.0f);
                float vawt_3d_correction = 1.0f / (1.0f + induction_factor);
                if (aspect_ratio < 2.0f) {
                    vawt_3d_correction *= (0.8f + 0.2f * aspect_ratio / 2.0f);
                }

                float3 f_blade = vawt_3d_correction * dyn_p_times_area * (ct * unit_tan + cn * unit_normal);
                float3 r_vec = pos - center;
                float torque_contribution = r_vec.x * f_blade.y - r_vec.y * f_blade.x;
                torque_lbm += torque_contribution;
                float3 f_fluid = -f_blade;

                // Distribute f_fluid to force field
                int range = (int)ceilf(6.0f * sigma);
                const float gauss_norm = 1.0f / pow(2.0f * M_PI * sigma * sigma, 1.5f);
                float total_weight = 0.0f;
                std::vector<std::tuple<int, int, int, float>> weighted_cells;
                for (int dx = -range; dx <= range; ++dx) {
                    for (int dy = -range; dy <= range; ++dy) {
                        for (int dz_g = -range; dz_g <= range; ++dz_g) {
                            int xx = (int)floorf(pos.x) + dx;
                            int yy = (int)floorf(pos.y) + dy;
                            int zz = (int)floorf(pos.z) + dz_g;
                            if (xx < 0 || xx >= (int)Nx || yy < 0 || yy >= (int)Ny || zz < 0 || zz >= (int)Nz) continue;
                            float3 grid_pos = float3((float)xx + 0.5f, (float)yy + 0.5f, (float)zz + 0.5f);
                            float3 dr = grid_pos - pos;
                            float r_squared = dot(dr, dr);
                            float g = gauss_norm * expf(-r_squared / (2.0f * sigma * sigma));
                            weighted_cells.push_back(std::make_tuple(xx, yy, zz, g));
                            total_weight += g;
                        }
                    }
                }

                if (total_weight > 1e-12f) {
                    float correction = 1.0f / total_weight;
                    for (const auto& [xx, yy, zz, g] : weighted_cells) {
                        ulong n = lbm.index((uint)xx, (uint)yy, (uint)zz);
                        float normalized_weight = g * correction;
                        lbm.F[n] += f_fluid.x * normalized_weight;
                        lbm.F[n + lbm.get_N()] += f_fluid.y * normalized_weight;
                        lbm.F[n + 2 * lbm.get_N()] += f_fluid.z * normalized_weight;
                    }
                }
            }
        }

        lbm.F.write_to_device();
        lbm.run(lbm_dt);
        theta += lbm_omega * (float)lbm_dt;
        if (theta > 2.0f * M_PI) theta -= 2.0f * M_PI;
        if (theta < -2.0f * M_PI) theta += 2.0f * M_PI;

        if (lbm.get_t() % 10 == 0 && lbm.get_t() >= 2000) {
            float si_torque = units.si_M(torque_lbm);
            float si_time = units.si_t((float)lbm.get_t());
            torque_file << si_time << "\t" << si_torque << "\n";
            torque_sum += si_torque;
            measurement_count++;
            torque_in_rev.push_back(si_torque);
            float delta_theta = theta - last_theta;
            if (delta_theta < 0) delta_theta += rev_angle;
            if (delta_theta >= rev_angle) {
                if (!torque_in_rev.empty()) {
                    float sum_rev = std::accumulate(torque_in_rev.begin(), torque_in_rev.end(), 0.0f);
                    float avg_rev = sum_rev / torque_in_rev.size();
                    rev_avg_torque_history.push_back(avg_rev);
                    if (rev_avg_torque_history.size() > history_size) {
                        rev_avg_torque_history.pop_front();
                    }
                    if (rev_avg_torque_history.size() == history_size) {
                        float sum = std::accumulate(rev_avg_torque_history.begin(), rev_avg_torque_history.end(), 0.0f);
                        float mean = sum / history_size;
                        float variance = 0.0f;
                        for (float avg_tq : rev_avg_torque_history) {
                            float diff = avg_tq - mean;
                            variance += diff * diff;
                        }
                        variance /= (history_size - 1);
                        float stddev = std::sqrt(variance);
                        if (mean != 0.0f && (stddev / std::abs(mean)) < convergence_threshold) {
                            converged = true;
                            printf("Convergence reached at time step %lu with coefficient of variation on rev averages = %f\n", lbm.get_t(), (stddev / std::abs(mean)));
                        }
                    }
                }
                torque_in_rev.clear();
                last_theta = theta;
            }

            lbm.u.read_from_device();
            auto measure_wake = [&](uint x_measure) -> float {
                float u_wake_sum = 0.0f;
                int valid_cells = 0;
                for (uint y = 0; y < Ny; y++) {
                    for (uint z = 0; z < Nz; z++) {
                        ulong n = lbm.index(x_measure, y, z);
                        float dy = (y - center.y) * dx;
                        float dz = (z - center.z) * dx;
                        float r = sqrtf(dy * dy + dz * dz);
                        if (r <= 2.0f * radius_si) {
                            float u_si_x = units.si_u(lbm.u[n]);
                            u_wake_sum += u_si_x;
                            valid_cells++;
                        }
                    }
                }
                float u_wake_avg = valid_cells > 0 ? u_wake_sum / valid_cells : 0.0f;
                if (u_wake_avg < 0.0f) u_wake_avg = 0.0f;
                return u_wake_avg;
            };

            float u_wake_3 = measure_wake(x_measure_3);
            float u_wake_5 = measure_wake(x_measure_5);
            float u_wake_10 = (x_measure_10 < Nx) ? measure_wake(x_measure_10) : -1.0f;
            float a = (si_velocity - u_wake_3) / (2.0f * si_velocity);
            if (a < 0.0f) a = 0.0f;
            if (a > 0.5f) a = 0.5f;
            float power_extracted = 4.0f * a * pow(1.0f - a, 2) * 0.5f * si_density * swept_area * pow(si_velocity, 3);
            float torque_wake = power_extracted / si_omega;
            float avg_torque = torque_sum / measurement_count;
            float avg_power = avg_torque * si_omega;

            float sum_u_si = 0.0f;
            ulong num_fluid_cells = 0ull;
            for (uint x = 0; x < (uint)center.x; x++) {
                for (uint y = 0; y < Ny; y++) {
                    for (uint z = 0; z < Nz; z++) {
                        ulong n = lbm.index(x, y, z);
                        if (lbm.flags[n] == 0u) {
                            sum_u_si += units.si_u(lbm.u[n]);
                            num_fluid_cells++;
                        }
                    }
                }
            }
            float avg_domain_u = (num_fluid_cells > 0) ? sum_u_si / num_fluid_cells : 0.0f;

            printf("Time %lu: ALM Torque = %.4f N·m, Avg ALM Torque = %.4f N·m, Avg Power = %.4f W, Wake Torque = %.4f N·m\n",
                   lbm.get_t(), si_torque, avg_torque, avg_power, torque_wake);
            printf("Wake Velocity at 3 radius (~1.35m downstream): %.4f m/s (expected <1.62 m/s)\n", u_wake_3);
            printf("Wake Velocity at 5 radius (~2.25m downstream): %.4f m/s (expected deficit, recovering)\n", u_wake_5);
            if (u_wake_10 >= 0.0f) {
                printf("Wake Velocity at 10 radius (~4.5m downstream): %.4f m/s (expected near 1.62 m/s)\n", u_wake_10);
            }
            printf("Average domain velocity (x-component, upstream only): %.4f m/s (expected ~1.62 m/s in steady state)\n", avg_domain_u);

            MLDataPoint data_point;
            data_point.time = si_time;
            data_point.theta = theta;
            data_point.torque = si_torque;
            data_point.TSR = (radius_si * fabsf(si_omega)) / si_velocity;
            data_point.Re_avg = si_velocity * chord_si / si_nu;
            data_point.inflow_velocity = avg_domain_u;
            data_point.wake_deficit = (si_velocity - u_wake_3) / si_velocity;
            ml_dataset.push_back(data_point);
        }

#if defined(GRAPHICS) && !defined(INTERACTIVE_GRAPHICS)
        if (lbm.graphics.next_frame(lbm_T, 30.0f)) {
            lbm.graphics.set_camera_free(float3(0.5f * Nx, 0.5f * Ny, 2.0f * Nz), 0.0f, 90.0f, 200.0f);
            lbm.graphics.write_frame();
        }
#endif
    }

    torque_file.close();
    std::ofstream ml_file("ml_data.csv");
    ml_file << "time,theta,torque,TSR,Re_avg,inflow_velocity,wake_deficit\n";
    for (const auto& dp : ml_dataset) {
        ml_file << dp.time << "," << dp.theta << "," << dp.torque << "," << dp.TSR << "," << dp.Re_avg << "," << dp.inflow_velocity << "," << dp.wake_deficit << "\n";
    }
    ml_file.close();

    auto end_time = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> elapsed = end_time - start_time;
    printf("Total simulation time in real world: %.2f seconds\n", elapsed.count());
}