diff options
Diffstat (limited to 'tangle')
66 files changed, 4189 insertions, 0 deletions
diff --git a/tangle/LLBM/base.h b/tangle/LLBM/base.h new file mode 100644 index 0000000..e356e02 --- /dev/null +++ b/tangle/LLBM/base.h @@ -0,0 +1,6 @@ +#pragma once + +#include "descriptor.h" +#include "memory.h" +#include "lattice.h" +#include "materials.h" diff --git a/tangle/LLBM/boundary.h b/tangle/LLBM/boundary.h new file mode 100644 index 0000000..08793a0 --- /dev/null +++ b/tangle/LLBM/boundary.h @@ -0,0 +1,13 @@ +#include "kernel/bounce_back.h" + +#include "kernel/bounce_back_moving_wall.h" + +#include "kernel/free_slip.h" + +#include "kernel/bouzidi.h" + +#include "kernel/equilibrium_velocity_wall.h" + +#include "kernel/equilibrium_density_wall.h" + +#include "sdf_boundary.h" diff --git a/tangle/LLBM/bulk.h b/tangle/LLBM/bulk.h new file mode 100644 index 0000000..62d704a --- /dev/null +++ b/tangle/LLBM/bulk.h @@ -0,0 +1,3 @@ +#include "kernel/collide.h" + +#include "kernel/smagorinsky_collide.h" diff --git a/tangle/LLBM/call_tag.h b/tangle/LLBM/call_tag.h new file mode 100644 index 0000000..fe559ef --- /dev/null +++ b/tangle/LLBM/call_tag.h @@ -0,0 +1,12 @@ +#pragma once + +namespace tag { + +struct call_by_cell_id { }; +struct call_by_list_index { }; +struct call_by_spatial_cell_mask { }; + +struct post_process_by_list_index { }; +struct post_process_by_spatial_cell_mask { }; + +} diff --git a/tangle/LLBM/descriptor.h b/tangle/LLBM/descriptor.h new file mode 100644 index 0000000..d4ca44d --- /dev/null +++ b/tangle/LLBM/descriptor.h @@ -0,0 +1,298 @@ +#pragma once + +#include <algorithm> +#include <cstdint> +#include <type_traits> +#include <cuda-samples/Common/helper_math.h> + +#ifdef __CUDA_ARCH__ + #define DATA device_data +#else + #define DATA host_data +#endif + +using pop_index_t = std::uint8_t; + +namespace descriptor { + +struct D2Q9 { + static constexpr unsigned d = 2; + static constexpr unsigned q = 9; +}; + +struct D3Q19 { + static constexpr unsigned d = 3; + static constexpr unsigned q = 19; +}; + +namespace device_data { + template <typename DESCRIPTOR> + __constant__ pop_index_t opposite[DESCRIPTOR::q] { }; + + template <typename DESCRIPTOR> + __constant__ int c[DESCRIPTOR::q][DESCRIPTOR::d] { }; + + template <typename DESCRIPTOR> + __constant__ float c_length[DESCRIPTOR::q] { }; + + template <typename DESCRIPTOR> + __constant__ float weight[DESCRIPTOR::q] { }; + + template <> + __constant__ pop_index_t opposite<D2Q9>[9] = { + 8, 7, 6, 5, 4, 3, 2, 1, 0 + }; + + template <> + __constant__ int c<D2Q9>[9][2] = { + {-1,-1}, {-1,0}, {-1,1}, {0,-1}, {0,0}, {0,1}, {1,-1}, {1,0}, {1,1} + }; + + template <> + __constant__ float c_length<D2Q9>[9] = { + 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.00000000000000, 0, 1.00000000000000, 1.41421356237310, 1.00000000000000, 1.41421356237310 + }; + + template <> + __constant__ float weight<D2Q9>[9] = { + 0.0277777777777778, 0.111111111111111, 0.0277777777777778, 0.111111111111111, 0.444444444444444, 0.111111111111111, 0.0277777777777778, 0.111111111111111, 0.0277777777777778 + }; + + template <> + __constant__ pop_index_t opposite<D3Q19>[19] = { + 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0 + }; + + template <> + __constant__ int c<D3Q19>[19][3] = { + {0,1,1}, {-1,0,1}, {0,0,1}, {1,0,1}, {0,-1,1}, {-1,1,0}, {0,1,0}, {1,1,0}, {-1,0,0}, {0,0,0}, {1,0,0}, {-1,-1,0}, {0,-1,0}, {1,-1,0}, {0,1,-1}, {-1,0,-1}, {0,0,-1}, {1,0,-1}, {0,-1,-1} + }; + + template <> + __constant__ float c_length<D3Q19>[19] = { + 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.00000000000000, 0, 1.00000000000000, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310 + }; + + template <> + __constant__ float weight<D3Q19>[19] = { + 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0555555555555556, 0.333333333333333, 0.0555555555555556, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778 + }; +} + +namespace host_data { + template <typename DESCRIPTOR> + constexpr pop_index_t opposite[DESCRIPTOR::q] { }; + + template <typename DESCRIPTOR> + constexpr int c[DESCRIPTOR::q][DESCRIPTOR::d] { }; + + template <typename DESCRIPTOR> + constexpr float c_length[DESCRIPTOR::q] { }; + + template <typename DESCRIPTOR> + constexpr float weight[DESCRIPTOR::q] { }; + + template <> + constexpr pop_index_t opposite<D2Q9>[9] = { + 8, 7, 6, 5, 4, 3, 2, 1, 0 + }; + + template <> + constexpr int c<D2Q9>[9][2] = { + {-1,-1}, {-1,0}, {-1,1}, {0,-1}, {0,0}, {0,1}, {1,-1}, {1,0}, {1,1} + }; + + template <> + constexpr float c_length<D2Q9>[9] = { + 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.00000000000000, 0, 1.00000000000000, 1.41421356237310, 1.00000000000000, 1.41421356237310 + }; + + template <> + constexpr float weight<D2Q9>[9] = { + 0.0277777777777778, 0.111111111111111, 0.0277777777777778, 0.111111111111111, 0.444444444444444, 0.111111111111111, 0.0277777777777778, 0.111111111111111, 0.0277777777777778 + }; + + template <> + constexpr pop_index_t opposite<D3Q19>[19] = { + 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0 + }; + + template <> + constexpr int c<D3Q19>[19][3] = { + {0,1,1}, {-1,0,1}, {0,0,1}, {1,0,1}, {0,-1,1}, {-1,1,0}, {0,1,0}, {1,1,0}, {-1,0,0}, {0,0,0}, {1,0,0}, {-1,-1,0}, {0,-1,0}, {1,-1,0}, {0,1,-1}, {-1,0,-1}, {0,0,-1}, {1,0,-1}, {0,-1,-1} + }; + + template <> + constexpr float c_length<D3Q19>[19] = { + 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.00000000000000, 0, 1.00000000000000, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310, 1.41421356237310, 1.00000000000000, 1.41421356237310, 1.41421356237310 + }; + + template <> + constexpr float weight<D3Q19>[19] = { + 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0555555555555556, 0.333333333333333, 0.0555555555555556, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778, 0.0277777777777778, 0.0555555555555556, 0.0277777777777778, 0.0277777777777778 + }; +} + +template <typename DESCRIPTOR> +__host__ __device__ +pop_index_t opposite(pop_index_t iPop) { + return DESCRIPTOR::q - 1 - iPop; +} + +template <typename DESCRIPTOR> +__host__ __device__ +int velocity(pop_index_t iPop, unsigned iDim) { + return DATA::template c<DESCRIPTOR>[iPop][iDim]; +} + +template <typename DESCRIPTOR> +__host__ __device__ +std::enable_if_t<DESCRIPTOR::d == 2, float2> velocity(pop_index_t iPop) { + return make_float2(DATA::template c<DESCRIPTOR>[iPop][0], + DATA::template c<DESCRIPTOR>[iPop][1]); +} + +template <typename DESCRIPTOR> +__host__ __device__ +std::enable_if_t<DESCRIPTOR::d == 3, float3> velocity(pop_index_t iPop) { + return make_float3(DATA::template c<DESCRIPTOR>[iPop][0], + DATA::template c<DESCRIPTOR>[iPop][1], + DATA::template c<DESCRIPTOR>[iPop][2]); +} + +template <typename DESCRIPTOR> +__host__ __device__ +float velocity_length(pop_index_t iPop) { + return DATA::template c_length<DESCRIPTOR>[iPop]; +} + +template <typename DESCRIPTOR> +__host__ __device__ +float weight(pop_index_t iPop) { + return DATA::template weight<DESCRIPTOR>[iPop]; +} + +template <unsigned D> +struct CuboidD; + +template <> +struct CuboidD<2> { + const std::size_t nX; + const std::size_t nY; + const std::size_t nZ; + const std::size_t volume; + + CuboidD(std::size_t x, std::size_t y): + nX(x), nY(y), nZ(1), + volume(x*y) { }; +}; + +template <> +struct CuboidD<3> { + const std::size_t nX; + const std::size_t nY; + const std::size_t nZ; + const std::size_t volume; + const std::size_t plane; + + CuboidD(std::size_t x, std::size_t y, std::size_t z): + nX(x), nY(y), nZ(z), + volume(x*y*z), + plane(x*y) { }; +}; + +template <typename DESCRIPTOR> +using Cuboid = CuboidD<DESCRIPTOR::d>; + +__host__ __device__ +std::size_t gid(const CuboidD<2>& c, int iX, int iY, int iZ=0) { + return iY*c.nX + iX; +} + +__host__ __device__ +std::size_t gid(const CuboidD<3>& c, int iX, int iY, int iZ) { + return iZ*c.plane + iY*c.nX + iX; +} + +__host__ __device__ +int offset(const CuboidD<2>& c, int iX, int iY) { + return iY*c.nX + iX; +} + +template <typename DESCRIPTOR> +__host__ __device__ +int offset(const CuboidD<2>& c, pop_index_t iPop) { + static_assert(DESCRIPTOR::d == 2, "Dimensions must match"); + return offset(c, + descriptor::velocity<DESCRIPTOR>(iPop, 0), + descriptor::velocity<DESCRIPTOR>(iPop, 1) + ); +} + +__host__ __device__ +int offset(const CuboidD<3>& c, int iX, int iY, int iZ) { + return iZ*c.plane + iY*c.nX + iX; +} + +template <typename DESCRIPTOR> +__host__ __device__ +int offset(const CuboidD<3>& c, pop_index_t iPop) { + static_assert(DESCRIPTOR::d == 3, "Dimensions must match"); + return offset(c, + descriptor::velocity<DESCRIPTOR>(iPop, 0), + descriptor::velocity<DESCRIPTOR>(iPop, 1), + descriptor::velocity<DESCRIPTOR>(iPop, 2) + ); +} + +template <typename DESCRIPTOR> +__host__ __device__ +std::size_t neighbor(const CuboidD<2>& c, std::size_t iCell, pop_index_t iPop) { + return iCell + offset<DESCRIPTOR>(c, iPop); +} + +template <typename DESCRIPTOR> +__host__ __device__ +std::size_t neighbor(const CuboidD<3>& c, std::size_t iCell, pop_index_t iPop) { + return iCell + offset<DESCRIPTOR>(c, iPop); +} + +__host__ __device__ +uint2 gidInverse(const CuboidD<2>& c, std::size_t gid) { + int iY = gid / c.nX; + int iX = gid % c.nX; + return make_uint2(iX, iY); +} + +__host__ __device__ +uint3 gidInverse(const CuboidD<3>& c, std::size_t gid) { + int iZ = gid / c.plane; + int iY = (gid % c.plane) / c.nX; + int iX = (gid % c.plane) % c.nX; + return make_uint3(iX,iY,iZ); +} + +__host__ __device__ +float2 gidInverseSmooth(const CuboidD<2>& c, std::size_t gid) { + int iY = gid / c.nX; + int iX = gid % c.nX; + return make_float2(iX, iY); +} + +__host__ __device__ +float3 gidInverseSmooth(const CuboidD<3>& c, std::size_t gid) { + int iZ = gid / c.plane; + int iY = (gid % c.plane) / c.nX; + int iX = (gid % c.plane) % c.nX; + return make_float3(iX,iY,iZ); +} + +bool isInside(const CuboidD<2>& c, std::size_t gid) { + return gid < c.volume; +} + +bool isInside(const CuboidD<3>& c, std::size_t gid) { + return gid < c.volume; +} + +} diff --git a/tangle/LLBM/kernel/bounce_back.h b/tangle/LLBM/kernel/bounce_back.h new file mode 100644 index 0000000..d7ec4a7 --- /dev/null +++ b/tangle/LLBM/kernel/bounce_back.h @@ -0,0 +1,44 @@ +#pragma once +#include <LLBM/call_tag.h> + +struct BounceBackO { + +using call_tag = tag::call_by_cell_id; + +template <typename T, typename S> +__device__ static void apply(descriptor::D2Q9, S f_curr[9], S f_next[9], std::size_t gid) { + f_next[0] = f_curr[8]; + f_next[1] = f_curr[7]; + f_next[2] = f_curr[6]; + f_next[3] = f_curr[5]; + f_next[4] = f_curr[4]; + f_next[5] = f_curr[3]; + f_next[6] = f_curr[2]; + f_next[7] = f_curr[1]; + f_next[8] = f_curr[0]; +} + +template <typename T, typename S> +__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid) { + f_next[0] = f_curr[18]; + f_next[1] = f_curr[17]; + f_next[2] = f_curr[16]; + f_next[3] = f_curr[15]; + f_next[4] = f_curr[14]; + f_next[5] = f_curr[13]; + f_next[6] = f_curr[12]; + f_next[7] = f_curr[11]; + f_next[8] = f_curr[10]; + f_next[9] = f_curr[9]; + f_next[10] = f_curr[8]; + f_next[11] = f_curr[7]; + f_next[12] = f_curr[6]; + f_next[13] = f_curr[5]; + f_next[14] = f_curr[4]; + f_next[15] = f_curr[3]; + f_next[16] = f_curr[2]; + f_next[17] = f_curr[1]; + f_next[18] = f_curr[0]; +} + +}; diff --git a/tangle/LLBM/kernel/bounce_back_moving_wall.h b/tangle/LLBM/kernel/bounce_back_moving_wall.h new file mode 100644 index 0000000..c4b40f9 --- /dev/null +++ b/tangle/LLBM/kernel/bounce_back_moving_wall.h @@ -0,0 +1,44 @@ +#pragma once +#include <LLBM/call_tag.h> + +struct BounceBackMovingWallO { + +using call_tag = tag::call_by_cell_id; + +template <typename T, typename S> +__device__ static void apply(descriptor::D2Q9, S f_curr[9], S f_next[9], std::size_t gid, T u_0, T u_1) { + f_next[0] = -T{0.166666666666667}*u_0 - T{0.166666666666667}*u_1 + f_curr[8]; + f_next[1] = -T{0.666666666666667}*u_0 + f_curr[7]; + f_next[2] = -T{0.166666666666667}*u_0 + T{0.166666666666667}*u_1 + f_curr[6]; + f_next[3] = -T{0.666666666666667}*u_1 + f_curr[5]; + f_next[4] = f_curr[4]; + f_next[5] = T{0.666666666666667}*u_1 + f_curr[3]; + f_next[6] = T{0.166666666666667}*u_0 - T{0.166666666666667}*u_1 + f_curr[2]; + f_next[7] = T{0.666666666666667}*u_0 + f_curr[1]; + f_next[8] = T{0.166666666666667}*u_0 + T{0.166666666666667}*u_1 + f_curr[0]; +} + +template <typename T, typename S> +__device__ static void apply(descriptor::D3Q19, S f_curr[19], S f_next[19], std::size_t gid, T u_0, T u_1, T u_2) { + f_next[0] = T{0.166666666666667}*u_1 + T{0.166666666666667}*u_2 + f_curr[18]; + f_next[1] = -T{0.166666666666667}*u_0 + T{0.166666666666667}*u_2 + f_curr[17]; + f_next[2] = T{0.333333333333333}*u_2 + f_curr[16]; + f_next[3] = T{0.166666666666667}*u_0 + T{0.166666666666667}*u_2 + f_curr[15]; + f_next[4] = -T{0.166666666666667}*u_1 + T{0.166666666666667}*u_2 + f_curr[14]; + f_next[5] = -T{0.166666666666667}*u_0 + T{0.166666666666667}*u_1 + f_curr[13]; + f_next[6] = T{0.333333333333333}*u_1 + f_curr[12]; + f_next[7] = T{0.166666666666667}*u_0 + T{0.166666666666667}*u_1 + f_curr[11]; + f_next[8] = -T{0.333333333333333}*u_0 + f_curr[10]; + f_next[9] = f_curr[9]; + f_next[10] = T{0.333333333333333}*u_0 + f_curr[8]; + f_next[11] = -T{0.166666666666667}*u_0 - T{0.166666666666667}*u_1 + f_curr[7]; + f_next[12] = -T{0.333333333333333}*u_1 + f_curr[6]; + f_next[13] = T{0.166666666666667}*u_0 - T{0.166666666666667}*u_1 + f_curr[5]; + f_next[14] = T{0.166666666666667}*u_1 - T{0.166666666666667}*u_2 + f_curr[4]; + f_next[15] = -T{0.166666666666667}*u_0 - T{0.166666666666667}*u_2 + f_curr[3]; + f_next[16] = -T{0.333333333333333}*u_2 + f_curr[2]; + f_next[17] = T{0.166666666666667}*u_0 - T{0.166666666666667}*u_2 + f_curr[1]; + f_next[18] = -T{0.166666666666667}*u_1 - T{0.166666666666667}*u_2 + f_curr[0]; +} + +}; diff --git a/tangle/LLBM/kernel/bouzidi.h b/tangle/LLBM/kernel/bouzidi.h new file mode 100644 index 0000000..b1abd2c --- /dev/null +++ b/tangle/LLBM/kernel/bouzidi.h @@ -0,0 +1,39 @@ +#pragma once +#include <LLBM/call_tag.h> +#include <LLBM/lattice.h> + +template <typename S> +struct BouzidiConfig { + std::size_t* boundary; // boundary cell to be interpolated + std::size_t* solid; // adjacent solid cell + std::size_t* fluid; // adjacent fluid cell + S* distance; // precomputed distance factor q + S* correction; // correction for moving walls + pop_index_t* missing; // population to be reconstructed +}; + +struct BouzidiO { + +using call_tag = tag::call_by_list_index; + +template <typename T, typename S, typename DESCRIPTOR> +__device__ static void apply( + LatticeView<DESCRIPTOR,S> lattice + , std::size_t index + , std::size_t count + , BouzidiConfig<S> config +) { + pop_index_t& iPop = config.missing[index]; + pop_index_t jPop = descriptor::opposite<DESCRIPTOR>(iPop); + pop_index_t kPop = config.boundary[index] == config.fluid[index] ? iPop : jPop; + + S f_bound_j = *lattice.pop(jPop, config.boundary[index]); + S f_fluid_j = *lattice.pop(kPop, config.fluid[index]); + S* f_next_i = lattice.pop(iPop, config.solid[index]); + + *f_next_i = config.distance[index] * f_bound_j + + (1. - config.distance[index]) * f_fluid_j + + config.correction[index]; +} + +}; diff --git a/tangle/LLBM/kernel/collect_curl.h b/tangle/LLBM/kernel/collect_curl.h new file mode 100644 index 0000000..054b359 --- /dev/null +++ b/tangle/LLBM/kernel/collect_curl.h @@ -0,0 +1,44 @@ +#pragma once +#include <LLBM/call_tag.h> + +struct CollectCurlF { + +using call_tag = tag::call_by_spatial_cell_mask; + +template <typename T, typename S> +__device__ static void apply( + descriptor::D3Q19 + , S f_curr[19] + , descriptor::CuboidD<3> cuboid + , std::size_t gid + , std::size_t iX + , std::size_t iY + , std::size_t iZ + , S* moments_u + , cudaSurfaceObject_t surface + , S* curl_norm = nullptr +) { + auto u_x = [moments_u,cuboid,gid] __device__ (int x, int y, int z) -> T { + return moments_u[3*(gid + descriptor::offset(cuboid,x,y,z)) + 0]; + }; + auto u_y = [moments_u,cuboid,gid] __device__ (int x, int y, int z) -> T { + return moments_u[3*(gid + descriptor::offset(cuboid,x,y,z)) + 1]; + }; + auto u_z = [moments_u,cuboid,gid] __device__ (int x, int y, int z) -> T { + return moments_u[3*(gid + descriptor::offset(cuboid,x,y,z)) + 2]; + }; + + T curl_0 = T{0.500000000000000}*u_y(0, 0, -1) - T{0.500000000000000}*u_y(0, 0, 1) - T{0.500000000000000}*u_z(0, -1, 0) + T{0.500000000000000}*u_z(0, 1, 0); + T curl_1 = -T{0.500000000000000}*u_x(0, 0, -1) + T{0.500000000000000}*u_x(0, 0, 1) + T{0.500000000000000}*u_z(-1, 0, 0) - T{0.500000000000000}*u_z(1, 0, 0); + T curl_2 = T{0.500000000000000}*u_x(0, -1, 0) - T{0.500000000000000}*u_x(0, 1, 0) - T{0.500000000000000}*u_y(-1, 0, 0) + T{0.500000000000000}*u_y(1, 0, 0); + float3 curl = make_float3(curl_0, curl_1, curl_2); + float norm = length(curl); + + surf3Dwrite(norm, surface, iX*sizeof(float), iY, iZ); + + if (curl_norm != nullptr) { + curl_norm[gid] = norm; + } +} + +}; diff --git a/tangle/LLBM/kernel/collect_moments.h b/tangle/LLBM/kernel/collect_moments.h new file mode 100644 index 0000000..b2ef09b --- /dev/null +++ b/tangle/LLBM/kernel/collect_moments.h @@ -0,0 +1,45 @@ +#pragma once +#include <LLBM/call_tag.h> + +struct CollectMomentsF { + +using call_tag = tag::call_by_cell_id; + +template <typename T, typename S> +__device__ static void apply(descriptor::D2Q9, S f_curr[9], std::size_t gid, T* cell_rho, T* cell_u) { + T m0 = f_curr[1] + f_curr[2]; + T m1 = f_curr[3] + f_curr[6]; + T m2 = m0 + m1 + f_curr[0] + f_curr[4] + f_curr[5] + f_curr[7] + f_curr[8]; + T m3 = f_curr[0] - f_curr[8]; + T m4 = T{1} / (m2); + T rho = m2; + T u_0 = -m4*(m0 + m3 - f_curr[6] - f_curr[7]); + T u_1 = -m4*(m1 + m3 - f_curr[2] - f_curr[5]); + + cell_rho[gid] = rho; + cell_u[2*gid+0] = u_0; + cell_u[2*gid+1] = u_1; +} + +template <typename T, typename S> +__device__ static void apply(descriptor::D3Q19, S f_curr[19], std::size_t gid, T* cell_rho, T* cell_u) { + T m0 = f_curr[10] + f_curr[13] + f_curr[17]; + T m1 = f_curr[14] + f_curr[5] + f_curr[6]; + T m2 = f_curr[1] + f_curr[2] + f_curr[4]; + T m3 = m0 + m1 + m2 + f_curr[0] + f_curr[11] + f_curr[12] + f_curr[15] + f_curr[16] + f_curr[18] + f_curr[3] + f_curr[7] + f_curr[8] + f_curr[9]; + T m4 = -f_curr[11] + f_curr[7]; + T m5 = -f_curr[15] + f_curr[3]; + T m6 = T{1} / (m3); + T m7 = f_curr[0] - f_curr[18]; + T rho = m3; + T u_0 = m6*(m0 + m4 + m5 - f_curr[1] - f_curr[5] - f_curr[8]); + T u_1 = m6*(m1 + m4 + m7 - f_curr[12] - f_curr[13] - f_curr[4]); + T u_2 = m6*(m2 + m5 + m7 - f_curr[14] - f_curr[16] - f_curr[17]); + + cell_rho[gid] = rho; + cell_u[3*gid+0] = u_0; + cell_u[3*gid+1] = u_1; + cell_u[3*gid+2] = u_2; +} + +}; diff --git a/tangle/LLBM/kernel/collect_q_criterion.h b/tangle/LLBM/kernel/collect_q_criterion.h new file mode 100644 index 0000000..19b7f68 --- /dev/null +++ b/tangle/LLBM/kernel/collect_q_criterion.h @@ -0,0 +1,104 @@ +#pragma once +#include <LLBM/call_tag.h> + +struct CollectQCriterionF { + +using call_tag = tag::call_by_spatial_cell_mask; + +template <typename T, typename S> +__device__ static void apply( + descriptor::D3Q19 + , S f_curr[19] + , descriptor::CuboidD<3> cuboid + , std::size_t gid + , std::size_t iX + , std::size_t iY + , std::size_t iZ + , T* cell_rho + , T* cell_u + , T* cell_curl_norm + , cudaSurfaceObject_t surface + , T* cell_q = nullptr +) { + const T rho = cell_rho[gid]; + const T u_0 = cell_u[3*gid + 0]; + const T u_1 = cell_u[3*gid + 1]; + const T u_2 = cell_u[3*gid + 2]; + + T x0 = T{72.0000000000000}*f_curr[5]; + T x1 = T{72.0000000000000}*f_curr[13]; + T x2 = T{6.00000000000000}*u_1; + T x3 = -x2; + T x4 = -u_0; + T x5 = x4 + u_1; + T x6 = u_1*u_1; + T x7 = T{3.00000000000000}*x6; + T x8 = u_0*u_0; |