// Actual production step kernels under a CPU launch shim. No CUDA claim.
#include <cmath>
#include <cstdio>
#include <cstdint>
#include <vector>
#define __device__
#define __constant__
#define __global__
#define __forceinline__ inline
#define __launch_bounds__(...)
#define LATTICE_D2Q9 0
#define LATTICE_D3Q19 1
#define LATTICE_MODEL 0
#define NQ 9
#define DIM 2
#define NZ 1
#define NT 1
#define COLLISION_MODEL 0
#define BLOCK_ROOT_BC 2
#define USE_LES 0
#define RHO 1.f
#define U0 .025f
#define OMEGA_COLLISION_MIN .01f
#define OMEGA_COLLISION_MAX 1.99f
#define Y_WALL_BC 0
#define OUTLET_MODE 0
struct Dim {unsigned int x=0,y=0,z=0;} threadIdx,blockIdx,blockDim;
inline float __fdiv_rn(float a,float b){return a/b;}
#include "core/precision.cuh"
#include "core/descriptors.cuh"
#include "core/layout.cuh"
#include "core/params.cuh"
#include "core/flags.cuh"
#include "operators/macro.cuh"
#include "operators/collision_srt.cuh"
#include "streaming/pull_double_buffer.cuh"
#if UNIFORM
#include "operators/forcing_guo.cuh"
#include "operators/helpers.cuh"
#include "boundary/bounce_back.cuh"
#include "boundary/inlet_outlet.cuh"
#include "step/one_step_double.cu"
#else
constexpr unsigned char KIND_OUTER=2;
#include "step/block_channel.cu"
#endif
int main(){
    constexpr int gx=24,gy=20,n=gx*gy;
    std::vector<float> state(9*n),next(9*n);
    for(int q=0;q<9;q++)for(int y=0;y<gy;y++)for(int x=0;x<gx;x++)
        state[q*n+y*gx+x]=d_w[q]*(1.f+.001f*std::sin(float(x+3*y+q)))-USE_DDF_SHIFTING*d_w[q];
    d_params={1.25f,0,0,0,0};
#if UNIFORM
    std::vector<uint16_t> flags(n,FLAG_FLUID);
    for(int y=0;y<gy;y++)for(int x=0;x<gx;x++){
        if(y==0||y==gy-1) flags[y*gx+x]=FLAG_SOLID|FLAG_BC_WALL;
        else if(x==0) flags[y*gx+x]=FLAG_SOLID|FLAG_BC_INLET;
        else if(x==gx-1) flags[y*gx+x]=FLAG_SOLID|FLAG_BC_OUTLET;
    }
#else
    constexpr int bx=NX-2,by=NY-2,nbx=gx/bx,nby=gy/by;
    std::vector<float> input(nbx*nby*9*TOTAL_CELLS),output(input.size());
    d_block_channel={gx,gy,U0};
#endif
    for(int step=0;step<5;step++){
#if UNIFORM
        for(int y=0;y<gy;y++)for(int x=0;x<gx;x++){
            blockIdx={unsigned(x),unsigned(y),0};OneStep(flags.data(),state.data(),next.data());
        }
#else
        // Populate cores and every wrapped halo from the same dense source.
        for(int b=0;b<nbx*nby;b++)for(int y=0;y<NY;y++)for(int x=0;x<NX;x++)for(int q=0;q<9;q++){
            int xx=((b%nbx)*bx+x-1+gx)%gx, yy=((b/nbx)*by+y-1+gy)%gy;
            input[b*9*TOTAL_CELLS+q*TOTAL_CELLS+y*NX+x]=state[q*n+yy*gx+xx];
        }
        for(int b=0;b<nbx*nby;b++)for(int y=1;y<NY-1;y++)for(int x=1;x<NX-1;x++){
            blockIdx={unsigned(x),unsigned(y),unsigned(b)};
            BlockChannelStep(input.data(),output.data(),1,nullptr,nullptr);
            for(int q=0;q<9;q++)next[q*n+((b/nbx)*by+y-1)*gx+(b%nbx)*bx+x-1]=output[b*9*TOTAL_CELLS+q*TOTAL_CELLS+y*NX+x];
        }
#endif
        state.swap(next);
    }
    for(float f:state)printf("%.9g\n",f);
}