#include <initializer_list>
#include <cstdlib>

#include <cmath>
#include <cstdio>
#include <cstring>
#include <cstdint>
#define __device__
#define __global__
#define __forceinline__ inline
#define __constant__
#define LATTICE_D2Q9 0
#define LATTICE_D3Q19 1
#define LATTICE_MODEL 0
#define DIM 2
#define NQ 9
#define NX 8
#define NY 8
#define NZ 1
#define NT 1
#define USE_LES 0
#define RHO 1.0f
#ifndef COLLISION_MODEL
#define COLLISION_MODEL 0
#endif
#define OMEGA_COLLISION_MIN .001f
#define OMEGA_COLLISION_MAX 1.99f
struct Idx {unsigned x=0,y=0,z=0;}threadIdx,blockIdx,blockDim;
struct Params{float omega=1.25f,fx=0,fy=0;}d_params;
inline float atomicAdd(float*p,float v){*p+=v;return *p;}
inline int obs_force_index(int,int){return 0;}
inline int obs_torque_index(int,int){return 0;}
#include "core/flags.cuh"
#include "core/descriptors.cuh"
#include "core/precision.cuh"
#include "core/layout.cuh"
#include "operators/macro.cuh"
#include "operators/collision_srt.cuh"
#include "operators/collision_trt.cuh"
#include "operators/collision_mrt.cuh"
#include "operators/forcing_guo.cuh"
#include "operators/helpers.cuh"
#include "streaming/esopull_single_buffer.cuh"
#include "streaming/esopull_semantic_helpers.cuh"
#include "boundary/bounce_back.cuh"
#include "boundary/curved_boundary.cuh"
#include "step/esopull_sync.cu"
#define U0 .03f
#define OBS_BODY_SLOT_FLOATS 7
#define __launch_bounds__(...)
inline int obs_sensor_index(int,int){return 0;}
#include "boundary/inlet_outlet.cuh"
#include "streaming/pull_double_buffer.cuh"
#include "step/one_step_double.cu"
#include "step/one_step_esopull.cu"
#include "step/aux_kernels.cu"
#include "step/init_flow.cu"
#include <cassert>
#include <algorithm>
void close(float a,float b,float tol=3e-7f) {
    if(fabsf(a-b)>tol) {fprintf(stderr,"got %.9g expected %.9g\n",a,b); abort();}
}
int main() {
    float fi[NX*NY*NQ],before[NX*NY*NQ],f[NQ],physical[NQ],staged[6*NX+2*NY*NQ];
    // Nonuniform populations: exact same-time address inverse, both parities.
    for(unsigned long t=0;t<4;t++) {
        for(int i=0;i<NX*NY*NQ;i++) fi[i]=.001f*((i*71)%113);
        memcpy(before,fi,sizeof fi);
        for(unsigned cell : {9u,27u,54u}) {
            unsigned long j[NQ]; compute_neighbors(cell,j);
            load_physical_node_esopull(cell,fi,j,t,f);
            store_physical_node_esopull(cell,f,fi,j,t);
            assert(memcmp(before,fi,sizeof fi)==0);
        }
        for(unsigned cell : {9u,27u,54u}) {
            // Decode/encode preserves arbitrary non-equilibrium, not just feq.
            DecodeCellsToPhysical(fi,physical,&cell,1,t);
            EncodePhysicalToCells(fi,physical,&cell,1,t);
        }
        for(int i=0;i<NX*NY*NQ;i++) close(fi[i],before[i],4e-8f);
    }
    // Stage y-wall donors from exact labelled postcollision states. Both
    // parities and wall rows. Staging is complete before any backing write.
    for(unsigned long t=1;t<3;t++) {
        for(unsigned long k=0;k<NX*NY;k++) {
            unsigned long j[NQ];compute_neighbors(k,j);
            for(int i=0;i<NQ;i++) f[i]=100*i+k;
            store_f_esopull(k,f,fi,j,t-1);
        }
        for(unsigned side=0;side<2;side++) for(unsigned x=0;x<NX;x++) {
            blockIdx.y=side;threadIdx.x=x;PrepareEsoPullBoundaries(fi,staged,t);
            unsigned k=(side?NY-2:1)*NX+x;
            unsigned donor[3]={side?3u:4u,side?5u:6u,side?8u:7u};
#if Y_WALL_BC == 1
            std::swap(donor[1],donor[2]);
#endif
            for(int a=0;a<3;a++) close(staged[(side*NX+x)*3+a],100*donor[a]+k,0);
        }
    }
    blockIdx.y=0;threadIdx.x=0;
    // t=0 uses the real initialization mapping, not double-buffer raw slots.
    for(unsigned long k=0;k<NX*NY;k++) {unsigned x,y;coordinates(k,x,y);write_equilibrium_esopull(fi,k,x,y,.03f);}
    compute_feq(1,.03f,0,f);
    for(unsigned side=0;side<2;side++) for(unsigned x=0;x<NX;x++) {
        blockIdx.y=side;threadIdx.x=x;PrepareEsoPullBoundaries(fi,staged,0);
        unsigned donor[3]={side?3u:4u,side?5u:6u,side?8u:7u};
#if Y_WALL_BC == 1
        std::swap(donor[1],donor[2]);
#endif
        for(int a=0;a<3;a++) close(staged[(side*NX+x)*3+a],f[donor[a]]);
    }
    blockIdx.y=0;threadIdx.x=0;

    // Complete open-boundary donor staging, including initial-state mapping,
    // odd/even backing, corners and both known-direction free-slip repairs.
    for(unsigned long t=0;t<3;t++) {
        for(unsigned long k=0;k<NX*NY;k++) {
            unsigned long j[NQ];compute_neighbors(k,j);
            for(unsigned i=0;i<NQ;i++)f[i]=.001f*(100*i+k);
            if(t)store_f_esopull(k,f,fi,j,t-1);
            else store_physical_node_esopull(k,f,fi,j,t);
        }
        for(unsigned side=0;side<2;side++) for(unsigned n=0;n<std::max(NX,NY);n++) {
            blockIdx.y=side;threadIdx.x=n;PrepareEsoPullBoundaries(fi,staged,t);
            if(n>=NY)continue;
            const float* donor=staged+6*NX+(side*NY+n)*NQ;
            unsigned long k=linear_index(side?NX-2:1,n);
            for(unsigned i=0;i<NQ;i++)close(donor[i],.001f*(100*i+k),0);
            float known[NQ]={0};repair_west_knowns_from_donor_d2q9(known,donor,0,n);
            if(n==1u){close(known[3],donor[3],0);close(known[8],donor[8],0);}
            if(n==NY-2u){close(known[4],donor[4],0);close(known[6],donor[6],0);}
        }
    }
    blockIdx.y=0;threadIdx.x=0;

#if USE_DDF_SHIFTING
    // Initialization must not form w*(1+tiny)-w and erase a representable flow.
    for(int layout=0;layout<2;layout++) {
        for(unsigned long k=0;k<NX*NY;k++) {
            const float velocity=1e-9f*(1+k%7);
            if(layout) {
                unsigned x,y;coordinates(k,x,y);
                write_equilibrium_esopull(fi,k,x,y,velocity);
            } else write_equilibrium(fi,k,velocity);
        }
        for(unsigned long k=0;k<NX*NY;k++) {
            unsigned long j[NQ];compute_neighbors(k,j);
            if(layout)load_f_esopull(k,f,fi,j,0);
            else for(unsigned i=0;i<NQ;i++)f[i]=fi[index_f(k,i)];
            float rho,u,v;compute_rho_u(f,rho,u,v);
            close(rho,1,1e-7f);close(u,1e-9f*(1+k%7),2e-15f);close(v,0,2e-15f);
        }
    }
#endif
    // Moving-wall equilibrium with rho != 1, all links/branches/fallbacks.
    for(float q : {.2f,.5f,.8f}) for(unsigned fallback : {0u,1u}) {
        compute_feq(1.17f,.02f,-.01f,f);
        for(unsigned dir=1;dir<NQ;dir++) {
            float value=compute_bouzidi_reflection(f[dir],f[dir],f[opp_dir(dir)],q,fallback)
                +compute_moving_wall_addend(dir,q,fallback,.02f,-.01f,1.17f);
            close(value,f[opp_dir(dir)]);
        }
    }

    // Verify both curved kernel call paths supply physical density to the
    // moving correction, not just the pure arithmetic helper above.
    for(unsigned long t=0;t<3;t++) {
        compute_feq(1.17f,.02f,-.01f,f);
        for(unsigned long k=0;k<NX*NY;k++) {
            unsigned long j[NQ];compute_neighbors(k,j);
            if(t)store_f_esopull(k,f,fi,j,t-1);
            else for(unsigned i=0;i<NQ;i++)fi[index_f(k,i)]=f[i];
        }
        unsigned cell=27;unsigned char dir=7,fallback=0,slot;
        int bid=0;float q=.5f,rx=-10,ry=-20,rz=0,corr;
        float act[OBS_BODY_SLOT_FLOATS]={0};act[OBS_BODY_SLOT_FLOATS-1]=.001f;
        if(t) {
            unsigned target;
            PrepareEsoPullCurvedKernel(fi,&cell,&dir,&q,&rx,&ry,&fallback,&bid,
                act,nullptr,&corr,&target,&slot,1,t);
        } else {
            CurvedBoundaryKernel(fi,&cell,&dir,&q,&rx,&ry,&rz,&fallback,&bid,act,nullptr,1);
            unsigned long j[NQ];compute_neighbors(cell,j);
            corr=fi[index_f(j[dir],opp_dir(dir))];
        }
        close(corr,f[opp_dir(dir)]);
    }
    // Execute actual region kernel: raw layout and same-time EsoPull, non-unit
    // density, nonzero u, opposite signs. Entire vector matches exact difference.
    unsigned cell=27;int bid=0;uint16_t flags[NX*NY];
    for(auto& fl:flags) fl=FLAG_FLUID;
    float action[7]={.001f,-.0007f,0,0,0,0,0};
    for(unsigned long t=0;t<3;t++) {
        float old[NQ],expected[NQ],new_eq[NQ];
        compute_feq(1.17f,.03f,-.02f,old);
        compute_feq(1.17f,.03f+action[0]/1.17f,-.02f+action[1]/1.17f,new_eq);
        for(int i=0;i<NQ;i++) expected[i]=new_eq[i];
        for(unsigned long k=0;k<NX*NY;k++) {
            unsigned long j[NQ];compute_neighbors(k,j);
            if(t)store_physical_node_esopull(k,old,fi,j,t);
            else for(int i=0;i<NQ;i++) fi[index_f(k,i)]=old[i];
        }
        ForceRegionKernel(fi,flags,&cell,&bid,action,1,t);
        unsigned long j[NQ];compute_neighbors(cell,j);
        if(t)load_physical_node_esopull(cell,fi,j,t,f);
        else for(int i=0;i<NQ;i++)f[i]=fi[index_f(cell,i)];
        for(int i=0;i<NQ;i++)close(f[i],expected[i]);
    }

    // Nonuniform closed-channel evolution through actual step kernels. x is
    // periodic, both y walls explicit; compare postcollision spatial values.
    float db[NX*NY*NQ],out[NX*NY*NQ];
    for(unsigned long k=0;k<NX*NY;k++) {
        unsigned x,y;coordinates(k,x,y);
        flags[k]=(y==0||y==NY-1)?FLAG_SOLID|FLAG_BC_WALL:FLAG_FLUID;
        compute_feq(1+.001f*cosf(k),.01f*sinf(x+y),.005f*cosf(x-y),f);
        unsigned long j[NQ];compute_neighbors(k,j);
        store_f_esopull(k,f,fi,j,0);
        for(int i=0;i<NQ;i++)db[index_f(k,i)]=f[i];
    }
    for(unsigned long t=1;t<=4;t++) {
        for(unsigned side=0;side<2;side++) for(unsigned x=0;x<NX;x++) {
            threadIdx.x=x;blockIdx.y=side;PrepareEsoPullBoundaries(fi,staged,t);
        }
        for(unsigned y=0;y<NY;y++) for(unsigned x=0;x<NX;x++) {
            threadIdx.x=x;blockIdx.y=y;
            OneStep(flags,db,out);
            EsoPullStep(fi,flags,nullptr,nullptr,staged,t);
        }
        for(unsigned y=1;y<NY-1;y++) for(unsigned x=0;x<NX;x++) {
            unsigned long k=linear_index(x,y),j[NQ];compute_neighbors(k,j);
            for(unsigned i=0;i<NQ;i++)
                close(load_physical_dir_pre_streaming(k,i,fi,j,t+1),out[index_f(k,i)]);
        }
        memcpy(db,out,sizeof db);
    }
    puts("same-time sync, staged wall donors, moving-wall and force-region passed");
}