#include <cmath>
#include <cstdio>
#include <initializer_list>
#define __device__
#define __constant__
#define __forceinline__ inline
#define LATTICE_D2Q9 0
#define LATTICE_MODEL 0
#define NQ 9
#define DIM 2
#define USE_DDF_SHIFTING 1
#define WALL_NORMAL_VELOCITY 0.f
#define WALL_SLOT_WIDTH 0
#define NX 16
#include "core/descriptors.cuh"
// Host round-to-nearest shim; CUDA intrinsic behavior is tested on-device.
inline float __fdiv_rn(float a, float b) { return a / b; }
#include "operators/macro.cuh"
#include "operators/collision_srt.cuh"
#include "boundary/planar_velocity.cuh"
void emit(const float* f) {for(int q=0;q<9;q++) printf("%.9g ",f[q]);puts("");}
int main() {
    for(float delta : {0.f, 0x1p-30f, -0x1p-30f, .0002f})
    for(float u : {0.f,.005f}) {
        float eq[9]; compute_feq(1.f+delta,u,-.2f*u,eq,delta);
        printf("E %.9g %.9g %.9g\n",delta,u,-.2f*u);emit(eq);
    }
    for(float delta : {0.f, 0x1p-30f, -0x1p-30f, .04f, -.01f})
    for(float jx : {0.f,.005f,-.02f}) {
        const float jy = -.2f*jx, rho = 1.f+delta;
        const float ux = jx/rho, uy = jy/rho;
        float eq[9];compute_feq_from_momentum(delta,jx,jy,ux,uy,eq);
        printf("J %.9g %.9g %.9g %.9g %.9g\n",delta,jx,jy,ux,uy);emit(eq);
    }
    for(int mode : {0,1,2}) for(float rho : {1.f,1.04f}) {
        float f[9];const float u=.007f,v=.003f;
        compute_feq(rho,u,v,f,rho-1.f);
        for(int q=0;q<9;q++)f[q]+=.00001f*(q-4);
        printf("B %d %.9g %.9g %.9g\n",mode,rho,u,v);emit(f);
        if(mode==0)reconstruct_bottom_velocity(f,u,v);
        else if(mode==1)reconstruct_top_velocity(f,u,v);
        else reconstruct_top_reservoir(f,rho);
        emit(f);
    }
}