aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--boltz.cc99
-rw-r--r--src/vector.h45
2 files changed, 72 insertions, 72 deletions
diff --git a/boltz.cc b/boltz.cc
index e4dd7b2..470d5d8 100644
--- a/boltz.cc
+++ b/boltz.cc
@@ -1,59 +1,22 @@
#include <iostream>
-#include <cmath>
-#include <algorithm>
#include <numeric>
#include <iostream>
#include <fstream>
-inline double sq(double x) noexcept {
- return x * x;
-}
+#include "src/vector.h"
struct DataCell {
double data[3][3];
- double& get(int x, int y) {
+ inline double& get(int x, int y) {
return data[1+x][1-y];
}
- double get(int x, int y) const {
+ inline double get(int x, int y) const {
return data[1+x][1-y];
}
};
-struct Vector {
- double data[2];
-
- double comp(int x, int y) const {
- return x*data[0] + y*data[1];
- }
-
- double norm() const {
- return std::sqrt(sq(data[0]) + sq(data[1]));
- }
-
- double& operator[](std::size_t i) {
- return data[i];
- }
-
- double operator[](std::size_t i) const {
- return data[i];
- }
-
- Vector operator*(double scalar) const {
- return Vector{
- data[0] * scalar,
- data[1] * scalar
- };
- }
-
- Vector& operator+=(const Vector& rhs) {
- data[0] += rhs[0];
- data[1] += rhs[1];
- return *this;
- }
-};
-
using Velocity = Vector;
using Force = Vector;
using Density = double;
@@ -66,15 +29,11 @@ constexpr DataCell weight{
struct Cell : DataCell {
void equilibrize(Density d, Velocity v) {
- get(-1, 1) = weight.get(-1, 1) * d * ( 1 + 3*v.comp(-1, 1) + 4.5*sq(v.comp(-1, 1)) - 1.5*sq(v.norm()) );
- get( 0, 1) = weight.get( 0, 1) * d * ( 1 + 3*v.comp( 0, 1) + 4.5*sq(v.comp( 0, 1)) - 1.5*sq(v.norm()) );
- get( 1, 1) = weight.get( 1, 1) * d * ( 1 + 3*v.comp( 1, 1) + 4.5*sq(v.comp( 1, 1)) - 1.5*sq(v.norm()) );
- get(-1, 0) = weight.get(-1, 0) * d * ( 1 + 3*v.comp(-1, 0) + 4.5*sq(v.comp(-1, 0)) - 1.5*sq(v.norm()) );
- get( 0, 0) = weight.get( 0, 0) * d * ( 1 - 1.5*sq(v.norm()) );
- get( 1, 0) = weight.get( 1, 0) * d * ( 1 + 3*v.comp( 1, 0) + 4.5*sq(v.comp( 1, 0)) - 1.5*sq(v.norm()) );
- get(-1,-1) = weight.get(-1,-1) * d * ( 1 + 3*v.comp(-1,-1) + 4.5*sq(v.comp(-1,-1)) - 1.5*sq(v.norm()) );
- get( 0,-1) = weight.get( 0,-1) * d * ( 1 + 3*v.comp( 0,-1) + 4.5*sq(v.comp( 0,-1)) - 1.5*sq(v.norm()) );
- get( 1,-1) = weight.get( 1,-1) * d * ( 1 + 3*v.comp( 1,-1) + 4.5*sq(v.comp( 1,-1)) - 1.5*sq(v.norm()) );
+ for ( int i = -1; i <= 1; ++i ) {
+ for ( int j = -1; j <= 1; ++j ) {
+ get(i,j) = weight.get(i,j) * d * (1 + 3*v.comp(i,j) + 4.5*sq(v.comp(i,j)) - 1.5*sq(v.norm()));
+ }
+ }
}
double sum() const {
@@ -82,18 +41,17 @@ struct Cell : DataCell {
}
Velocity velocity(Density d) const {
- return Velocity{
+ return 1./d * Velocity{
get( 1, 0) - get(-1, 0) + get( 1, 1) - get(-1,-1) + get( 1,-1) - get(-1,1),
get( 0, 1) - get( 0,-1) + get( 1, 1) - get(-1,-1) - get( 1,-1) + get(-1,1)
- } * ( 1./d );
+ };
}
};
-constexpr std::size_t dimX = 128;
-constexpr std::size_t dimY = 128;
+constexpr std::size_t dimX = 128;
+constexpr std::size_t dimY = 128;
-constexpr double tau = 0.6;
-constexpr double omega = 1. / tau;
+constexpr double tau = 1.0;
Cell cellsA[dimX][dimY];
Cell cellsB[dimX][dimY];
@@ -121,28 +79,22 @@ 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 ) {
// update velocity (force coupling)
- //velocity[x][y] += force[x][y] * (tau / density[x][y]);
+ //velocity[x][y] += tau / density[x][y] * force[x][y];
// compute equilibrium
Cell eq;
eq.equilibrize(density[x][y], velocity[x][y]);
- Cell& old = (*oldCells)[x][y];
-
// collide & stream
- (*newCells)[x ][y ].get( 0, 0) = old.get( 0, 0) + omega * ( eq.get( 0, 0) - old.get( 0, 0) );
- (*newCells)[x + 1][y ].get( 1, 0) = old.get( 1, 0) + omega * ( eq.get( 1, 0) - old.get( 1, 0) );
- (*newCells)[x - 1][y ].get(-1, 0) = old.get(-1, 0) + omega * ( eq.get(-1, 0) - old.get(-1, 0) );
- (*newCells)[x ][y + 1].get( 0, 1) = old.get( 0, 1) + omega * ( eq.get( 0, 1) - old.get( 0, 1) );
- (*newCells)[x ][y - 1].get( 0,-1) = old.get( 0,-1) + omega * ( eq.get( 0,-1) - old.get( 0,-1) );
- (*newCells)[x + 1][y + 1].get( 1, 1) = old.get( 1, 1) + omega * ( eq.get( 1, 1) - old.get( 1, 1) );
- (*newCells)[x - 1][y - 1].get(-1,-1) = old.get(-1,-1) + omega * ( eq.get(-1,-1) - old.get(-1,-1) );
- (*newCells)[x + 1][y - 1].get( 1,-1) = old.get( 1,-1) + omega * ( eq.get( 1,-1) - old.get( 1,-1) );
- (*newCells)[x - 1][y + 1].get(-1, 1) = old.get(-1, 1) + omega * ( eq.get(-1, 1) - old.get(-1, 1) );
+ for ( int i = -1; i <= 1; ++i ) {
+ for ( int j = -1; j <= 1; ++j ) {
+ (*newCells)[x+i][y+j].get(i,j) = (*oldCells)[x][y].get(i,j) + 1./tau * (eq.get(i,j) - (*oldCells)[x][y].get(i,j));
+ }
+ }
}
}
- // straight wall cell bounce back
+ // straight wall cell bounce back
for ( std::size_t x = 2; x < dimX - 2; ++x ) {
(*newCells)[x][1].get( 0, 1) = (*newCells)[x][0].get( 0,-1);
(*newCells)[x][1].get( 1, 1) = (*newCells)[x][0].get(-1,-1);
@@ -242,8 +194,11 @@ void init() {
}
}
- for ( std::size_t x = 50; x < dimX-50; ++x ) {
- for ( std::size_t y = 50; y < dimY-50; ++y ) {
+ for ( std::size_t y = 55; y < dimY-55; ++y ) {
+ for ( std::size_t x = 35; x < dimX-75; ++x ) {
+ density[x][y] = 0.8;
+ }
+ for ( std::size_t x = 75; x < dimX-35; ++x ) {
density[x][y] = 0.8;
}
}
@@ -252,10 +207,10 @@ void init() {
int main() {
init();
- for ( std::size_t t = 0; t < 500; ++t ) {
+ for ( std::size_t t = 0; t < 100; ++t ) {
computeLbmStep(t);
- if ( t % 1 == 0 ) {
+ if ( t % 2 == 0 ) {
writeCurrentStateAsVTK(t);
}
}
diff --git a/src/vector.h b/src/vector.h
new file mode 100644
index 0000000..55c7663
--- /dev/null
+++ b/src/vector.h
@@ -0,0 +1,45 @@
+#include <cmath>
+
+inline double sq(double x) noexcept {
+ return x * x;
+}
+
+struct Vector {
+ double data[2];
+
+ double comp(int x, int y) const {
+ return x*data[0] + y*data[1];
+ }
+
+ double norm() const {
+ return std::sqrt(sq(data[0]) + sq(data[1]));
+ }
+
+ double& operator[](std::size_t i) {
+ return data[i];
+ }
+
+ double operator[](std::size_t i) const {
+ return data[i];
+ }
+
+ Vector operator*(double scalar) const {
+ return Vector{
+ data[0] * scalar,
+ data[1] * scalar
+ };
+ }
+
+ Vector& operator+=(const Vector& rhs) {
+ data[0] += rhs[0];
+ data[1] += rhs[1];
+ return *this;
+ }
+};
+
+Vector operator*(double scalar, const Vector& vec) {
+ return Vector{
+ vec[0] * scalar,
+ vec[1] * scalar
+ };
+}