94 lines
3.3 KiB
C++
94 lines
3.3 KiB
C++
// 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);
|
|
}
|