aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--boltz.cc113
-rw-r--r--src/vector.h22
2 files changed, 78 insertions, 57 deletions
diff --git a/boltz.cc b/boltz.cc
index 85d939a..1acebc2 100644
--- a/boltz.cc
+++ b/boltz.cc
@@ -18,8 +18,8 @@ struct DataCell {
}
};
-using Velocity = Vector;
-using Force = Vector;
+using Velocity = Vector<double>;
+using Force = Vector<double>;
using Density = double;
constexpr DataCell weight{
@@ -81,7 +81,7 @@ class CellBuffer {
constexpr std::size_t dimX = 128;
constexpr std::size_t dimY = 128;
-constexpr double tau = 0.6;
+constexpr double tau = 0.8;
CellBuffer pop(dimX, dimY);
@@ -90,6 +90,8 @@ Velocity velocity[dimX][dimY];
Force force [dimX][dimY];
void computeFluidCell(std::size_t x, std::size_t y) {
+ //velocity[x][y] += tau / density[x][y] * force[x][y];
+
// compute equilibrium
Cell eq;
eq.equilibrize(density[x][y], velocity[x][y]);
@@ -102,14 +104,31 @@ void computeFluidCell(std::size_t x, std::size_t y) {
}
}
-inline int clamp(int x) {
- return x / -x;
+std::pair<Vector<int>, Vector<int>> neighbors(Vector<int> v) {
+ if ( v[0] == 0 ) {
+ return std::pair<Vector<int>, Vector<int>>{
+ { -1, v[1] },
+ { 1, v[1] }
+ };
+ } else if ( v[1] == 0 ) {
+ return std::pair<Vector<int>, Vector<int>>{
+ { v[0], -1 },
+ { v[0], 1 }
+ };
+ } else {
+ return std::pair<Vector<int>, Vector<int>>{
+ { 0, v[1] },
+ { v[0], 0 }
+ };
+ }
}
-void computeWallCell(std::size_t x, std::size_t y, int normalX, int normalY) {
- pop.curr(x,y).get(clamp(normalX+normalY),clamp(normalY-normalX)) = pop.curr(x-normalX,y-normalY).get(clamp(normalX-normalY),clamp(normalY+normalX));
- pop.curr(x,y).get(normalX ,normalY ) = pop.curr(x-normalX,y-normalY).get(-normalX,-normalY);
- pop.curr(x,y).get(clamp(normalX-normalY),clamp(normalY+normalX)) = pop.curr(x-normalX,y-normalY).get(clamp(normalX+normalY),clamp(normalY-normalX));
+void computeWallCell(std::size_t x, std::size_t y, Vector<int> normal) {
+ const auto [neighborA, neighborB] = neighbors(normal);
+
+ pop.curr(x,y).get(neighborA[0], neighborA[1]) = pop.curr(x-neighborA[0], y-neighborA[1]).get(-neighborA[0], -neighborA[1]);
+ pop.curr(x,y).get(normal[0] , normal[1] ) = pop.curr(x-normal[0] , y-normal[1] ).get(-normal[0] , -normal[1] );
+ pop.curr(x,y).get(neighborB[0], neighborB[1]) = pop.curr(x-neighborB[0], y-neighborB[1]).get(-neighborB[0], -neighborB[1]);
}
void computeLbmStep(std::size_t t) {
@@ -117,7 +136,7 @@ void computeLbmStep(std::size_t t) {
for ( std::size_t x = 1; x < dimX - 1; ++x ) {
for ( std::size_t y = 1; y < dimY - 1; ++y ) {
- if ( x <= 20 || x >= 40 || y <= 20 || y >= 40 ) {
+ if ( x <= 20 || x >= 40 || y <= 20 || y >= 80 ) {
computeFluidCell(x, y);
}
}
@@ -126,35 +145,35 @@ void computeLbmStep(std::size_t t) {
// obstacle
{
for ( std::size_t x = 21; x < 40; ++x ) {
- computeWallCell(x, 40, 0, 1);
- computeWallCell(x, 20, 0, -1);
+ computeWallCell(x, 80, { 0, 1});
+ computeWallCell(x, 20, { 0,-1});
}
- for ( std::size_t y = 21; y < 40; ++y ) {
- computeWallCell(40, y, 1, 0);
- computeWallCell(20, y, -1, 0);
+ for ( std::size_t y = 21; y < 80; ++y ) {
+ computeWallCell(40, y, { 1, 0});
+ computeWallCell(20, y, {-1, 0});
}
- computeWallCell(40,40, 1, 1);
- computeWallCell(40,20, 1,-1);
- computeWallCell(20,40,-1, 1);
- computeWallCell(20,20,-1,-1);
+ computeWallCell(40, 80, { 1, 1});
+ computeWallCell(40, 20, { 1,-1});
+ computeWallCell(20, 80, {-1, 1});
+ computeWallCell(20, 20, {-1,-1});
}
// straight wall cell bounce back
for ( std::size_t x = 2; x < dimX - 2; ++x ) {
- computeWallCell(x, 1, 0, 1);
- computeWallCell(x, dimY-2, 0, -1);
+ computeWallCell(x, 1, { 0, 1});
+ computeWallCell(x, dimY-2,{ 0,-1});
}
for ( std::size_t y = 2; y < dimY - 2; ++y ) {
- computeWallCell(1,y,1,0);
- computeWallCell(dimX-2,y,-1,0);
+ computeWallCell(1, y, { 1, 0});
+ computeWallCell(dimX-2, y, {-1, 0});
}
// edge wall cell bounce back
- computeWallCell(1,1,1,1);
- computeWallCell(1,dimY-2,1,-1);
- computeWallCell(dimX-2,1,-1,1);
- computeWallCell(dimX-2,dimY-2,-1,-1);
+ computeWallCell(1, 1, { 1, 1});
+ computeWallCell(1, dimY-2, { 1,-1});
+ computeWallCell(dimX-2, 1, {-1, 1});
+ computeWallCell(dimX-2, dimY-2, {-1,-1});
// update density, velocity field
for ( std::size_t x = 0; x < dimX; ++x ) {
@@ -165,6 +184,25 @@ void computeLbmStep(std::size_t t) {
}
}
+void init() {
+ for ( std::size_t x = 0; x < dimX; ++x ) {
+ for ( std::size_t y = 0; y < dimY; ++y ) {
+ density[x][y] = 1.0;
+ velocity[x][y] = { 0.0, 0.0 };
+ force[x][y] = { 0.0, 0.0 };
+
+ pop.curr(x,y).equilibrize(density[x][y], velocity[x][y]);
+ pop.prev(x,y).equilibrize(density[x][y], velocity[x][y]);
+ }
+ }
+
+ for ( std::size_t y = 50; y < dimY-50; ++y ) {
+ for ( std::size_t x = 50; x < dimX-50; ++x ) {
+ density[x][y] = 2.0;
+ }
+ }
+}
+
void writeCurrentStateAsVTK(int time) {
std::ofstream fout;
fout.open(("result/data_t" + std::to_string(time) + ".vtk").c_str());
@@ -207,29 +245,10 @@ void writeCurrentStateAsVTK(int time) {
fout.close();
}
-void init() {
- for ( std::size_t x = 0; x < dimX; ++x ) {
- for ( std::size_t y = 0; y < dimY; ++y ) {
- density[x][y] = 1.0;
- velocity[x][y] = { 0.0, 0.0 };
- force[x][y] = { 0.0, 0.0 };
-
- pop.curr(x,y).equilibrize(density[x][y], velocity[x][y]);
- pop.prev(x,y).equilibrize(density[x][y], velocity[x][y]);
- }
- }
-
- for ( std::size_t y = 55; y < dimY-55; ++y ) {
- for ( std::size_t x = 75; x < dimX-35; ++x ) {
- density[x][y] = 0.6;
- }
- }
-}
-
int main() {
init();
- for ( std::size_t t = 0; t < 500; ++t ) {
+ for ( std::size_t t = 0; t < 1000; ++t ) {
computeLbmStep(t);
if ( t % 2 == 0 ) {
diff --git a/src/vector.h b/src/vector.h
index 55c7663..e0dc239 100644
--- a/src/vector.h
+++ b/src/vector.h
@@ -4,41 +4,43 @@ inline double sq(double x) noexcept {
return x * x;
}
+template <typename T>
struct Vector {
- double data[2];
+ T data[2];
- double comp(int x, int y) const {
+ T comp(int x, int y) const {
return x*data[0] + y*data[1];
}
- double norm() const {
+ T norm() const {
return std::sqrt(sq(data[0]) + sq(data[1]));
}
- double& operator[](std::size_t i) {
+ T& operator[](std::size_t i) {
return data[i];
}
- double operator[](std::size_t i) const {
+ T operator[](std::size_t i) const {
return data[i];
}
- Vector operator*(double scalar) const {
- return Vector{
+ Vector<T> operator*(T scalar) const {
+ return Vector<T>{
data[0] * scalar,
data[1] * scalar
};
}
- Vector& operator+=(const Vector& rhs) {
+ Vector<T>& operator+=(const Vector<T>& rhs) {
data[0] += rhs[0];
data[1] += rhs[1];
return *this;
}
};
-Vector operator*(double scalar, const Vector& vec) {
- return Vector{
+template <typename T>
+Vector<T> operator*(T scalar, const Vector<T>& vec) {
+ return Vector<T>{
vec[0] * scalar,
vec[1] * scalar
};