aboutsummaryrefslogtreecommitdiff
path: root/cavity2d.cc
blob: 0006252ce6160b51d2c2f5317bc8ef4335d3b07a (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
#include <iostream>
#include <vector>
#include <string>
#include <algorithm>

#include "lbm.h"
#include "boundary_conditions.h"

constexpr std::size_t dimX = 64;
constexpr std::size_t dimY = dimX;

constexpr double uLid     = 0.1;
constexpr double reynolds = 100;

constexpr double tau   = 3 * uLid * (dimX-1) / reynolds + 0.5;
constexpr double omega = 1. / tau;

DataCellBuffer pop(dimX, dimY);
FluidBuffer fluid(dimX, dimY);

void init() {
	for ( std::size_t x = 0; x < dimX; ++x ) {
		for ( std::size_t y = 0; y < dimY; ++y ) {
			fluid.density(x,y)  = 1.0;
			fluid.velocity(x,y) = { 0.0, 0.0 };

			static_cast<Cell&>(pop.curr(x,y)).equilibrize(
				fluid.density(x,y), fluid.velocity(x,y));
			static_cast<Cell&>(pop.prev(x,y)).equilibrize(
				fluid.density(x,y), fluid.velocity(x,y));
		}
	}
}

void computeLbmStep() {
	pop.swap();

	for ( std::size_t x = 0; x < dimX; ++x ) {
		for ( std::size_t y = 0; y < dimY; ++y ) {
			streamFluidCell(pop, x, y);
		}
	}

	// straight wall cell bounce back
	for ( std::size_t x = 0; x < dimX; ++x ) {
		computeZouHeVelocityWallCell(pop, {x, dimY-1}, { 0,-1}, uLid);
	}
	for ( std::size_t x = 1; x < dimX-1; ++x ) {
		computeWallCell(pop, {x, 0}, { 0, 1});
	}
	for ( std::size_t y = 1; y < dimY-1; ++y ) {
		computeWallCell(pop, {0,      y}, { 1, 0});
		computeWallCell(pop, {dimX-1, y}, {-1, 0});
	}

	// edge wall cell bounce back
	computeWallCell(pop, {0,      0 }, { 1, 1});
	computeWallCell(pop, {dimX-1, 0 }, {-1, 1});

	for ( std::size_t x = 0; x < dimX; ++x ) {
		for ( std::size_t y = 0; y < dimY; ++y ) {
			Cell& cell = static_cast<Cell&>(pop.curr(x,y));
			fluid.density(x,y)  = cell.sum();
			fluid.velocity(x,y) = cell.velocity(fluid.density(x,y));

			collideFluidCell(omega, pop, fluid, x, y);
		}
	}
}

int main() {
	init();

	std::cout << "Re:   " << reynolds << std::endl;
	std::cout << "uLid: " << uLid     << std::endl;
	std::cout << "tau:  " << tau      << std::endl;

	for ( std::size_t t = 0; t <= 6000; ++t ) {
		computeLbmStep();

		if ( t % 100 == 0 ) {
			std::cout << ".";
			std::cout.flush();
			fluid.writeAsVTK("result/data_t" + std::to_string(t) + ".vtk");
		}
	}

	std::cout << std::endl;
}