// CPU qualifier shim: executes actual operators, never creates a GPU context.
#include <cmath>
#include <initializer_list>
#include <cstdio>
#define __device__
#define __constant__
#define __forceinline__ inline
#define LATTICE_D2Q9 0
#define LATTICE_D3Q19 1
#define NQ (DIM == 2 ? 9 : 19)
#define USE_LES 0
#define OMEGA_COLLISION_MIN .001f
#define OMEGA_COLLISION_MAX 1.99f
struct Params {float omega, fx, fy, fz;} d_params;
#include "core/descriptors.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 "boundary/inlet/common.cuh"
void emit(const float* f) {for(int i=0;i<NQ;++i) printf("%.9g ",f[i]); puts("");}
int main() {
#ifdef MRT_WEAK_FORCE_CHECK
    // Exercise actual FP32 collision against an independent FP64 matrix oracle.
    for (float omega : {.7f, 1.25f, 1.8f})
    for (float density : {1.f, 1.000001f, 1.04f})
    for (float velocity : {0.f, .0075f})
    for (float force : {0.f, 1.46484375e-6f, 1.e-8f}) {
        d_params = {omega, force, -.3f * force, 0.f};
        float f[9], source[9], rho, u, v;
        compute_feq(density, velocity, -.4f * velocity, f);
        if (velocity != 0.f)
            for (int i=0; i<9; ++i) f[i] += 1.e-5f * ((i*7)%5-2);
        compute_rho_u(f, rho, u, v);
        float uc=u, vc=v;
        apply_guo_velocity_correction(uc, vc, d_params.fx, d_params.fy, rho);
        compute_guo_forcing(uc, vc, d_params.fx, d_params.fy, source);
        printf("%.9g %.9g %.9g %.9g %.9g %.9g\n", omega, rho, uc, vc,
               d_params.fx, d_params.fy);
        emit(f); emit(source);
        collide_dispatch(f, rho, u, v);
        emit(f);
    }
    return 0;
#endif
    for (float omega : {0.7f, 1.25f, 1.8f}) {
        d_params = {omega, .001f, -.0007f, .0003f};
        float f[NQ], feq[NQ], source[NQ];
#if DIM == 2
        compute_feq(1.17f,.03f,-.02f,f);
#else
        compute_feq(1.17f,.03f,-.02f,.01f,f);
#endif
        // Non-equilibrium content: the expected result uses actual raw moments.
        for(int i=0;i<NQ;++i) f[i] += .0001f * ((i*7)%5-2);
        float rho,u,v,w=0;
#if DIM == 2
        compute_rho_u(f,rho,u,v);
        float uc=u,vc=v;
        apply_guo_velocity_correction(uc,vc,d_params.fx,d_params.fy,rho);
        compute_feq(rho,uc,vc,feq);
        compute_guo_forcing(uc,vc,d_params.fx,d_params.fy,source);
#else
        compute_rho_u(f,rho,u,v,w);
        float uc=u,vc=v,wc=w;
        apply_guo_velocity_correction(uc,vc,wc,d_params.fx,d_params.fy,d_params.fz,rho);
        compute_feq(rho,uc,vc,wc,feq);
        compute_guo_forcing(uc,vc,wc,d_params.fx,d_params.fy,d_params.fz,source);
#endif
        emit(f); emit(feq); emit(source);
#if DIM == 2
        collide_dispatch(f,rho,u,v);
#else
        collide_dispatch(f,rho,u,v,w);
#endif
        emit(f);
    }
    // Equilibrium and shifted inlet closure for non-unit density.
    d_params={1.25f,0,0,0};
    float f[NQ];
#if DIM == 2
    compute_feq(1.17f,.03f,-.02f,f);
    printf("closure %.9g\n", west_velocity_rho_closure_d2q9(f,.03f));
    collide_dispatch(f,1.17f,.03f,-.02f);
#else
    compute_feq(1.17f,.03f,-.02f,.01f,f);
    printf("closure %.9g\n", west_velocity_rho_closure_d3q19(f,.03f));
    collide_dispatch(f,1.17f,.03f,-.02f,.01f);
#endif
    emit(f);
}