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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
/*
* Lattice Boltzmann grid refinement sample, written in C++,
* using the OpenLB library
*
* Copyright (C) 2019 Adrian Kummerländer
* E-mail contact: info@openlb.net
* The most recent release of OpenLB can be downloaded at
* <http://www.openlb.net/>
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public
* License along with this program; if not, write to the Free
* Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*/
#include "olb2D.h"
#include "olb2D.hh"
using namespace olb;
typedef double T;
#define DESCRIPTOR descriptors::D2Q9<>
/// Setup geometry relative to cylinder diameter as defined by [SchaeferTurek96]
const T cylinderD = 0.1;
const int N = 20; // resolution of the cylinder
const T Re = 100.; // Reynolds number
const T tau = 0.51; // relaxation time
const T maxPhysT = 16.; // max. simulation time in s, SI unit
const Characteristics<T> PhysCharacteristics(
cylinderD, // char. phys. length
1.0, // char. phys. velocity
0.1/Re, // phsy. kinematic viscosity
1.0); // char. phys. density
#include "../common/model.h"
void setupRefinement(Grid2D<T,DESCRIPTOR>& coarseGrid,
Vector<T,2> domainOrigin, Vector<T,2> domainExtend,
Vector<T,2> cylinderCenter)
{
const auto coarseDeltaX = coarseGrid.getConverter().getPhysDeltaX();
const Vector<T,2> fineExtend {11.2*cylinderD, domainExtend[1]-1.5*coarseDeltaX};
const Vector<T,2> fineOrigin {0.5*cylinderD, coarseDeltaX};
auto& fineGrid = coarseGrid.refine(fineOrigin, fineExtend);
SchaeferTurek::prepareGeometry(fineGrid, domainOrigin, domainExtend);
const Vector<T,2> fineExtend2 {5*cylinderD, fineGrid.getExtend()[1]-2*coarseDeltaX};
const Vector<T,2> fineOrigin2 {1*cylinderD, (domainExtend[1]-fineExtend2[1])/2};
auto& fineGrid2 = fineGrid.refine(fineOrigin2, fineExtend2);
SchaeferTurek::prepareGeometry(fineGrid2, domainOrigin, domainExtend);
const Vector<T,2> fineExtend3 {1.4*cylinderD, 1.4*cylinderD};
const Vector<T,2> fineOrigin3 {cylinderCenter[0]-fineExtend3[0]/2, cylinderCenter[1]-fineExtend3[1]/2};
auto& fineGrid3 = fineGrid2.refine(fineOrigin3, fineExtend3);
SchaeferTurek::prepareGeometry(fineGrid3, domainOrigin, domainExtend);
}
int main(int argc, char* argv[])
{
olbInit(&argc, &argv);
singleton::directories().setOutputDir("./tmp/");
OstreamManager clout(std::cout,"main");
IndicatorCuboid2D<T> coarseCuboid(SchaeferTurek::modelExtend, SchaeferTurek::modelOrigin);
Grid2D<T,DESCRIPTOR> coarseGrid(
coarseCuboid,
RelaxationTime<T>(tau),
N,
PhysCharacteristics);
const Vector<T,2> domainOrigin = coarseGrid.getSuperGeometry().getStatistics().getMinPhysR(0);
const Vector<T,2> domainExtend = coarseGrid.getSuperGeometry().getStatistics().getPhysExtend(0);
SchaeferTurek::prepareGeometry(coarseGrid, domainOrigin, domainExtend);
setupRefinement(coarseGrid, domainOrigin, domainExtend, SchaeferTurek::cylinderCenter);
coarseGrid.forEachGrid(SchaeferTurek::prepareLattice);
clout << "Total number of active cells: " << coarseGrid.getActiveVoxelN() << endl;
clout << "Starting simulation..." << endl;
const int statIter = coarseGrid.getConverter().getLatticeTime(0.01);
Timer<T> timer(
coarseGrid.getConverter().getLatticeTime(maxPhysT),
coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
timer.start();
Grid2D<T,DESCRIPTOR>& cylinderGrid = coarseGrid.locate(SchaeferTurek::cylinderCenter);
for (int iT = 0; iT <= coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
SchaeferTurek::setBoundaryValues(coarseGrid, iT);
coarseGrid.collideAndStream();
if (iT == 0 || iT%statIter == 0) {
timer.update(iT);
timer.printStep();
coarseGrid.forEachGrid("cylinder2d", [&](Grid2D<T,DESCRIPTOR>& grid, const std::string& id) {
SchaeferTurek::getResults(grid, id, iT);
});
SchaeferTurek::takeMeasurements(cylinderGrid, iT);
}
}
timer.stop();
timer.printSummary();
}
|
