Line data Source code
1 : // Copyright (C) 2015 Technische Universitaet Muenchen
2 : // This file is part of the Mamico project. For conditions of distribution
3 : // and use, please see the copyright notice in Mamico's main folder
4 : #ifndef _MOLECULARDYNAMICS_COUPLING_SOLVERS_LBCOUETTESOLVER_H_
5 : #define _MOLECULARDYNAMICS_COUPLING_SOLVERS_LBCOUETTESOLVER_H_
6 :
7 : namespace coupling {
8 : namespace solvers {
9 : class LBCouetteSolver;
10 : class LBCouetteSolverState;
11 : } // namespace solvers
12 : } // namespace coupling
13 :
14 : #if defined(_OPENMP)
15 : #include <omp.h>
16 : #endif
17 : #include "coupling/interface/PintableMacroSolver.h"
18 : #include "coupling/solvers/NumericalSolver.h"
19 : #include <cmath>
20 :
21 64 : class coupling::solvers::LBCouetteSolverState : public coupling::interface::PintableMacroSolverState {
22 : public:
23 5756 : LBCouetteSolverState(int size) : _pdf(size, 0) {}
24 :
25 98 : LBCouetteSolverState(int size, double* pdf) : LBCouetteSolverState(size) { std::copy(pdf, pdf + size, _pdf.data()); }
26 :
27 4 : std::unique_ptr<State> clone() const override { return std::make_unique<LBCouetteSolverState>(*this); }
28 :
29 192 : ~LBCouetteSolverState() {}
30 :
31 4 : int getSizeBytes() const override { return sizeof(double) * _pdf.size(); }
32 :
33 : std::unique_ptr<State> operator+(const State& rhs) override;
34 : std::unique_ptr<State> operator-(const State& rhs) override;
35 :
36 52 : bool operator==(const State& rhs) const override {
37 52 : const LBCouetteSolverState* other = dynamic_cast<const LBCouetteSolverState*>(&rhs);
38 52 : if (other == nullptr)
39 : return false;
40 52 : return _pdf == other->_pdf;
41 : }
42 :
43 53320 : double* getData() override { return _pdf.data(); }
44 16 : const double* getData() const override { return _pdf.data(); }
45 :
46 0 : void print(std::ostream& os) const override { os << "<LBCouetteSolverState instance with size " << getSizeBytes() << ">"; }
47 :
48 : private:
49 : std::vector<double> _pdf;
50 : };
51 :
52 : /** In our scenario, the lower wall is accelerated and the upper wall stands
53 : * still. The lower wall is located at zero height.
54 : * @brief implements a three-dimensional Lattice-Boltzmann Couette flow solver.
55 : * @author Philipp Neumann */
56 : class coupling::solvers::LBCouetteSolver : public coupling::solvers::NumericalSolver, public coupling::interface::PintableMacroSolver {
57 : public:
58 : /** @brief a simple constructor
59 : * @param channelheight the width and height of the channel in y and z
60 : * direction
61 : * @param wallVelocity velocity at the moving wall, refers to Couette
62 : * scenario
63 : * @param dx the spacial step size, and equidistant grid is applied
64 : * @param dt the time step
65 : * @param kinVisc the kinematic viscosity of the fluid
66 : * @param plotEveryTimestep the time step interval for plotting data;
67 : * 4 means, every 4th time step is plotted
68 : * @param plotAverageVelocity writes average velocity for all time steps into CSV file
69 : * @param filestem the name of the plotted file
70 : * @param processes defines on how many processes the solver will run;
71 : * 1,1,1 - sequential run - 1,2,2 = 1*2*2 = 4 processes
72 : * @param numThreads number of OpenMP threads */
73 116 : LBCouetteSolver(const double channelheight, tarch::la::Vector<3, double> wallVelocity, const double kinVisc, const double dx, const double dt,
74 : const int plotEveryTimestep, const bool plotAverageVelocity, const std::string filestem, const tarch::la::Vector<3, unsigned int> processes,
75 : const unsigned int numThreads = 1, const Scenario* scen = nullptr)
76 116 : : coupling::solvers::NumericalSolver(channelheight, dx, dt, kinVisc, plotEveryTimestep, filestem, processes, scen), _mode(Mode::coupling), _dt_pint(dt),
77 232 : _omega(1.0 / (3.0 * (kinVisc * dt / (dx * dx)) + 0.5)), _wallVelocity((dt / dx) * wallVelocity), _plotAverageVelocity(plotAverageVelocity) {
78 : // return if required
79 116 : if (skipRank()) {
80 : return;
81 : }
82 29 : _pdfsize = 19 * (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2);
83 29 : _pdf1 = new double[_pdfsize];
84 29 : _pdf2 = new double[_pdfsize];
85 : #if defined(_OPENMP)
86 29 : omp_set_num_threads(numThreads);
87 : #endif
88 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
89 : std::cout << "Domain size=" << _domainSizeX << "," << _domainSizeY << "," << _domainSizeZ << std::endl;
90 : std::cout << "tau=" << 1.0 / _omega << std::endl;
91 : std::cout << "wallVelocity=" << _wallVelocity << std::endl;
92 : for (int z = 0; z < _domainSizeZ + 2; z++) {
93 : for (int y = 0; y < _domainSizeY + 2; y++) {
94 : for (int x = 0; x < _domainSizeX + 2; x++) {
95 : std::cout << x << "," << y << "," << z << "FLAG=" << _flag[get(x, y, z)] << std::endl;
96 : }
97 : }
98 : }
99 : #endif
100 : // check pointers
101 29 : if ((_pdf1 == NULL) || (_pdf2 == NULL) || (_vel == NULL) || (_density == NULL) || (_flag == NULL)) {
102 0 : std::cout << "ERROR LBCouetteSolver: NULL ptr!" << std::endl;
103 0 : exit(EXIT_FAILURE);
104 : }
105 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
106 29 : if ((_sendBufferX == NULL) || (_recvBufferX == NULL) || (_sendBufferY == NULL) || (_recvBufferY == NULL) || (_sendBufferZ == NULL) ||
107 29 : (_recvBufferZ == NULL)) {
108 0 : std::cout << "ERROR LBCouetteSolver: NULL ptr in send/recv!" << std::endl;
109 0 : exit(EXIT_FAILURE);
110 : }
111 : #endif
112 : // init everything with lattice weights
113 29 : #pragma omp parallel for
114 : for (int i = 0; i < (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2); i++) {
115 : for (int q = 0; q < 19; q++) {
116 : _pdf1[get(i) * 19 + q] = _W[q];
117 : _pdf2[get(i) * 19 + q] = _W[q];
118 : }
119 : }
120 29 : computeDensityAndVelocityEverywhere();
121 0 : }
122 :
123 : /** @brief a simple destructor */
124 232 : virtual ~LBCouetteSolver() {
125 116 : if (_pdf1 != NULL) {
126 29 : delete[] _pdf1;
127 29 : _pdf1 = NULL;
128 : }
129 116 : if (_pdf2 != NULL) {
130 29 : delete[] _pdf2;
131 29 : _pdf2 = NULL;
132 : }
133 116 : if (_vel != NULL) {
134 116 : delete[] _vel;
135 116 : _vel = NULL;
136 : }
137 116 : if (_density != NULL) {
138 116 : delete[] _density;
139 116 : _density = NULL;
140 : }
141 116 : if (_flag != NULL) {
142 116 : delete[] _flag;
143 116 : _flag = NULL;
144 : }
145 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
146 116 : if (_sendBufferX != NULL) {
147 116 : delete[] _sendBufferX;
148 116 : _sendBufferX = NULL;
149 : }
150 116 : if (_sendBufferY != NULL) {
151 116 : delete[] _sendBufferY;
152 116 : _sendBufferY = NULL;
153 : }
154 116 : if (_sendBufferZ != NULL) {
155 116 : delete[] _sendBufferZ;
156 116 : _sendBufferZ = NULL;
157 : }
158 116 : if (_recvBufferX != NULL) {
159 116 : delete[] _recvBufferX;
160 116 : _recvBufferX = NULL;
161 : }
162 116 : if (_recvBufferY != NULL) {
163 116 : delete[] _recvBufferY;
164 116 : _recvBufferY = NULL;
165 : }
166 116 : if (_recvBufferZ != NULL) {
167 116 : delete[] _recvBufferZ;
168 116 : _recvBufferZ = NULL;
169 : }
170 : #endif
171 232 : }
172 :
173 : /** @brief advances one time step dt in time and triggers vtk plot if required
174 : */
175 36 : void advance(double dt) override {
176 36 : if (skipRank()) {
177 : return;
178 : }
179 12 : const int timesteps = floor(dt / _dt + 0.5);
180 12 : if (fabs(timesteps * _dt - dt) / _dt > 1.0e-8) {
181 0 : std::cout << "ERROR LBCouetteSolver::advance(): time steps and dt do not match!" << std::endl;
182 0 : exit(EXIT_FAILURE);
183 : }
184 63 : for (int i = 0; i < timesteps; i++) {
185 51 : if (_plotEveryTimestep >= 1 && _counter % _plotEveryTimestep == 0)
186 0 : computeDensityAndVelocityEverywhere();
187 51 : plot();
188 51 : plot_avg_vel();
189 51 : collidestream();
190 51 : communicate(); // exchange between neighbouring MPI subdomains
191 51 : _counter++;
192 : }
193 : }
194 :
195 : /** @brief applies the values received from the MD-solver within the
196 : * conntinuum solver
197 : * @param md2macroBuffer holds the data from the md solver
198 : * coupling cells */
199 0 : void setMDBoundaryValues(coupling::datastructures::FlexibleCellContainer<3>& md2macroBuffer) override {
200 0 : if (skipRank()) {
201 : return;
202 : }
203 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
204 0 : if (_mode == Mode::supervising) {
205 0 : std::cout << "ERROR LBCouetteSolver setMDBoundaryValues() called in supervising mode" << std::endl;
206 0 : exit(EXIT_FAILURE);
207 : }
208 : #endif
209 0 : computeDensityAndVelocityEverywhere();
210 :
211 : // loop over all received cells
212 0 : for (auto pair : md2macroBuffer) {
213 0 : I01 idx;
214 0 : const coupling::datastructures::CouplingCell<3>* couplingCell;
215 0 : std::tie(couplingCell, idx) = pair;
216 : // determine cell index of this cell in LB domain
217 0 : tarch::la::Vector<3, unsigned int> globalCellCoords{idx.get()};
218 0 : globalCellCoords[0] = (globalCellCoords[0] + _offset[0]) - _coords[0] * _avgDomainSizeX;
219 0 : globalCellCoords[1] = (globalCellCoords[1] + _offset[1]) - _coords[1] * _avgDomainSizeY;
220 0 : globalCellCoords[2] = (globalCellCoords[2] + _offset[2]) - _coords[2] * _avgDomainSizeZ;
221 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
222 : std::cout << "Process coords: " << _coords << ": GlobalCellCoords for index " << idx << ": " << globalCellCoords << std::endl;
223 : #endif
224 0 : const int index = get(globalCellCoords[0], globalCellCoords[1], globalCellCoords[2]);
225 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
226 0 : if (_flag[index] != MD_BOUNDARY) {
227 0 : std::cout << "ERROR LBCouetteSolver::setMDBoundaryValues(): Cell " << index << " is no MD boundary cell!" << std::endl;
228 0 : exit(EXIT_FAILURE);
229 : }
230 : #endif
231 : // set velocity value and pdfs in MD boundary cell (before streaming); the
232 : // boundary velocities are interpolated between the neighbouring and this
233 : // cell. This interpolation is valid for FLUID-MD_BOUNDARY neighbouring
234 : // relations only. determine local velocity received from MaMiCo and
235 : // convert it to LB units; store the velocity in _vel
236 : // massFactor is used to ensure conservation of energy
237 0 : tarch::la::Vector<3, double> localVel((1.0 / couplingCell->getMacroscopicMass()) * (_dt / _dx) * couplingCell->getMacroscopicMomentum());
238 0 : for (unsigned int d = 0; d < 3; d++) {
239 0 : _vel[3 * index + d] = localVel[d];
240 : }
241 : // loop over all pdfs and set them according to interpolated moving-wall
242 : // conditions
243 0 : for (unsigned int q = 0; q < 19; q++) {
244 : // index of neighbour cell; only if cell is located inside local domain
245 0 : if (((int)globalCellCoords[0] + _C[q][0] > 0) && ((int)globalCellCoords[0] + _C[q][0] < _domainSizeX + 1) &&
246 0 : ((int)globalCellCoords[1] + _C[q][1] > 0) && ((int)globalCellCoords[1] + _C[q][1] < _domainSizeY + 1) &&
247 0 : ((int)globalCellCoords[2] + _C[q][2] > 0) && ((int)globalCellCoords[2] + _C[q][2] < _domainSizeZ + 1)) {
248 0 : const int nbIndex = get((_C[q][0] + globalCellCoords[0]), (_C[q][1] + globalCellCoords[1]), (_C[q][2] + globalCellCoords[2]));
249 0 : const tarch::la::Vector<3, double> interpolVel(0.5 * (_vel[3 * index] + _vel[3 * nbIndex]), 0.5 * (_vel[3 * index + 1] + _vel[3 * nbIndex + 1]),
250 0 : 0.5 * (_vel[3 * index + 2] + _vel[3 * nbIndex + 2]));
251 0 : _pdf1[19 * index + q] =
252 0 : _pdf1[19 * nbIndex + 18 - q] -
253 0 : 6.0 * _W[q] * _density[nbIndex] * (_C[18 - q][0] * interpolVel[0] + _C[18 - q][1] * interpolVel[1] + _C[18 - q][2] * interpolVel[2]);
254 : }
255 : }
256 : }
257 : }
258 :
259 : /** @brief returns velocity at a certain position
260 : * @param pos position for which the velocity will be returned
261 : * @returns the velocity vector for the position */
262 300 : tarch::la::Vector<3, double> getVelocity(tarch::la::Vector<3, double> pos) const override {
263 300 : tarch::la::Vector<3, unsigned int> coords;
264 300 : const tarch::la::Vector<3, double> domainOffset(_coords[0] * _dx * _avgDomainSizeX, _coords[1] * _dx * _avgDomainSizeY, _coords[2] * _dx * _avgDomainSizeZ);
265 : // check pos-data for process locality (todo: put this in debug mode in
266 : // future releases)
267 300 : if ((pos[0] < domainOffset[0]) || (pos[0] > domainOffset[0] + _domainSizeX * _dx) || (pos[1] < domainOffset[1]) ||
268 600 : (pos[1] > domainOffset[1] + _domainSizeY * _dx) || (pos[2] < domainOffset[2]) || (pos[2] > domainOffset[2] + _domainSizeZ * _dx)) {
269 0 : std::cout << "ERROR LBCouetteSolver::getVelocity(): Position " << pos << " out of range!" << std::endl;
270 0 : std::cout << "domainOffset = " << domainOffset << std::endl;
271 0 : std::cout << "_domainSizeX = " << _domainSizeX << std::endl;
272 0 : std::cout << "_domainSizeY = " << _domainSizeY << std::endl;
273 0 : std::cout << "_domainSizeZ = " << _domainSizeZ << std::endl;
274 0 : std::cout << "_dx = " << _dx << std::endl;
275 0 : exit(EXIT_FAILURE);
276 : }
277 : // compute index for respective cell (_dx+... for ghost cells); use coords
278 : // to store local cell coordinates
279 1200 : for (unsigned int d = 0; d < 3; d++) {
280 900 : coords[d] = (unsigned int)((_dx + pos[d] - domainOffset[d]) / _dx);
281 : }
282 300 : const int index = get(coords[0], coords[1], coords[2]);
283 300 : tarch::la::Vector<3, double> vel(0.0);
284 : // extract and scale velocity to "real"=MD units
285 1200 : for (int d = 0; d < 3; d++) {
286 900 : vel[d] = _dx / _dt * _vel[3 * index + d];
287 : }
288 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
289 : std::cout << "Position " << pos << " corresponds to cell: " << coords << "; vel=" << vel << std::endl;
290 : #endif
291 300 : return vel;
292 : }
293 :
294 : /** @brief returns density at a certain position
295 : * @param pos position for which the density will be returned
296 : * @returns the density vector for the position */
297 300 : double getDensity(tarch::la::Vector<3, double> pos) const override {
298 300 : tarch::la::Vector<3, unsigned int> coords;
299 300 : const tarch::la::Vector<3, double> domainOffset(_coords[0] * _dx * _avgDomainSizeX, _coords[1] * _dx * _avgDomainSizeY, _coords[2] * _dx * _avgDomainSizeZ);
300 : // check pos-data for process locality (todo: put this in debug mode in
301 : // future releases)
302 300 : if ((pos[0] < domainOffset[0]) || (pos[0] > domainOffset[0] + _domainSizeX * _dx) || (pos[1] < domainOffset[1]) ||
303 600 : (pos[1] > domainOffset[1] + _domainSizeY * _dx) || (pos[2] < domainOffset[2]) || (pos[2] > domainOffset[2] + _domainSizeZ * _dx)) {
304 0 : std::cout << "ERROR LBCouetteSolver::getDensity(): Position " << pos << " out of range!" << std::endl;
305 0 : exit(EXIT_FAILURE);
306 : }
307 : // compute index for respective cell (_dx+... for ghost cells); use coords
308 : // to store local cell coordinates
309 1200 : for (unsigned int d = 0; d < 3; d++) {
310 900 : coords[d] = (unsigned int)((_dx + pos[d] - domainOffset[d]) / _dx);
311 : }
312 300 : const int index = get(coords[0], coords[1], coords[2]);
313 300 : return _density[index];
314 : }
315 :
316 : /** @brief changes the velocity at the moving wall (z=0)
317 : * @param wallVelocity the velocity will be set at the moving wall */
318 0 : virtual void setWallVelocity(const tarch::la::Vector<3, double> wallVelocity) override { _wallVelocity = (_dt / _dx) * wallVelocity; }
319 :
320 : /// ------------------------------------------------------------------------------------------------------
321 : /// Pint methods ------ Pint methods ------ Pint methods ------ Pint methods ------ Pint methods
322 : /// ------------------------------------------------------------------------------------------------------
323 :
324 80 : std::unique_ptr<State> getState() override {
325 80 : computeDensityAndVelocityEverywhere();
326 80 : if (skipRank())
327 54 : return std::make_unique<LBCouetteSolverState>(0);
328 26 : return std::make_unique<LBCouetteSolverState>(_pdfsize, _pdf1);
329 : }
330 :
331 23 : void setState(const std::unique_ptr<State>& input, int cycle) override {
332 23 : if (skipRank())
333 : return;
334 :
335 8 : const LBCouetteSolverState* state = dynamic_cast<const LBCouetteSolverState*>(input.get());
336 :
337 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
338 8 : if (state == nullptr) {
339 0 : std::cout << "ERROR LBCouetteSolver setState() wrong state type" << std::endl;
340 0 : exit(EXIT_FAILURE);
341 : }
342 : #endif
343 :
344 8 : std::copy(state->getData(), state->getData() + _pdfsize, _pdf1);
345 8 : computeDensityAndVelocityEverywhere();
346 :
347 8 : _counter = cycle;
348 : }
349 :
350 14 : std::unique_ptr<State> operator()(const std::unique_ptr<State>& input, int cycle) override {
351 14 : setState(input, cycle);
352 :
353 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
354 14 : if (_mode != Mode::supervising) {
355 0 : std::cout << "ERROR LBCouetteSolver operator() called but not in supervising mode" << std::endl;
356 0 : exit(EXIT_FAILURE);
357 : }
358 : #endif
359 :
360 14 : advance(_dt_pint);
361 14 : return getState();
362 : }
363 :
364 12 : Mode getMode() const override { return _mode; }
365 :
366 : /**
367 : * This will create a new instance of this LBCouetteSolver
368 : * In the supervisor, setMDBoundary() has not been called, independently of state of couette solver in coupling mode
369 : * The supervisor can run with a modified viscosity
370 : * The supervisor's operator() advances many coupling cycles at once, interval is stored in _dt_pint
371 : */
372 72 : std::unique_ptr<PintableMacroSolver> getSupervisor(int num_cycles, double visc_multiplier) const override {
373 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
374 72 : if (_mode == Mode::supervising) {
375 0 : std::cout << "ERROR LBCouetteSolver getSupervisor(): already in supervising mode" << std::endl;
376 0 : exit(EXIT_FAILURE);
377 : }
378 : #endif
379 :
380 72 : int numThreads = 1;
381 : #if defined(_OPENMP)
382 72 : numThreads = omp_get_num_threads();
383 : #endif
384 :
385 144 : auto res = std::make_unique<LBCouetteSolver>(_channelheight, _wallVelocity * _dx / _dt, _kinVisc * visc_multiplier, _dx, _dt, _plotEveryTimestep,
386 216 : _plotAverageVelocity, _filestem + std::string("_supervising"), _processes, numThreads, _scen);
387 :
388 72 : res->_mode = Mode::supervising;
389 72 : res->_dt_pint = _dt * num_cycles;
390 :
391 144 : return res;
392 72 : }
393 :
394 0 : void print(std::ostream& os) const override {
395 0 : if (_mode == Mode::supervising)
396 0 : os << "<LBCouetteSolver instance in supervising mode >";
397 0 : if (_mode == Mode::coupling)
398 0 : os << "<LBCouetteSolver instance in coupling mode >";
399 0 : }
400 :
401 29 : double get_avg_vel(const std::unique_ptr<State>& state) const override {
402 29 : if (skipRank())
403 : return 0;
404 5 : double vel[3];
405 5 : double density;
406 5 : double res[3]{0, 0, 0};
407 53245 : for (int i = 0; i < _pdfsize; i += 19) {
408 53240 : LBCouetteSolver::computeDensityAndVelocity(vel, density, state->getData() + i);
409 53240 : res[0] += vel[0];
410 53240 : res[1] += vel[1];
411 53240 : res[2] += vel[2];
412 : }
413 5 : if (_pdfsize > 0) {
414 5 : res[0] /= (_pdfsize / 19);
415 5 : res[1] /= (_pdfsize / 19);
416 5 : res[2] /= (_pdfsize / 19);
417 : }
418 5 : return std::sqrt(res[0] * res[0] + res[1] * res[1] + res[2] * res[2]);
419 : }
420 :
421 0 : double get_avg_velX(const std::unique_ptr<State>& state) const {
422 0 : if (skipRank())
423 : return 0;
424 : double vel[3];
425 : double density;
426 : double res{0};
427 0 : for (int i = 0; i < _pdfsize; i += 19) {
428 0 : LBCouetteSolver::computeDensityAndVelocity(vel, density, state->getData() + i);
429 0 : res += vel[0];
430 : }
431 0 : if (_pdfsize > 0) {
432 0 : res /= (_pdfsize / 19);
433 : }
434 : return res;
435 : }
436 :
437 : private:
438 : Mode _mode;
439 : double _dt_pint;
440 :
441 117 : void computeDensityAndVelocityEverywhere() {
442 117 : if (skipRank())
443 : return;
444 1323 : for (int z = 1; z < _domainSizeZ + 1; z++) {
445 26460 : for (int y = 1; y < _domainSizeY + 1; y++) {
446 529200 : for (int x = 1; x < _domainSizeX + 1; x++) {
447 504000 : const int index = get(x, y, z);
448 504000 : const int pI = 19 * index;
449 504000 : double* vel = &_vel[3 * index];
450 504000 : computeDensityAndVelocity(vel, _density[index], &_pdf1[pI]);
451 : }
452 : }
453 : }
454 : }
455 :
456 51 : void plot_avg_vel() {
457 51 : if (!_plotAverageVelocity)
458 51 : return;
459 :
460 0 : int rank = 0;
461 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
462 0 : MPI_Comm_rank(coupling::indexing::IndexingService<3>::getInstance().getComm(), &rank);
463 : #endif
464 0 : std::stringstream ss;
465 0 : ss << _filestem << "_r" << rank;
466 0 : if (_scen != nullptr) {
467 0 : auto ts = _scen->getTimeIntegrationService();
468 0 : if (ts != nullptr) {
469 0 : if (ts->isPintEnabled())
470 0 : ss << "_i" << ts->getIteration();
471 : }
472 : }
473 0 : ss << ".csv";
474 0 : std::string filename = ss.str();
475 0 : std::ofstream file(filename.c_str(), _counter == 0 ? std::ofstream::out : std::ofstream::app);
476 0 : if (!file.is_open()) {
477 0 : std::cout << "ERROR LBCouetteSolver::plot_avg_vel(): Could not open file " << filename << "!" << std::endl;
478 0 : exit(EXIT_FAILURE);
479 : }
480 :
481 0 : if (_counter == 0) {
482 0 : file << "coupling_cycle ; avg_vel ; avg_velX" << std::endl;
483 : }
484 :
485 0 : std::unique_ptr<State> s = std::make_unique<LBCouetteSolverState>(_pdfsize, _pdf1);
486 0 : double vel = get_avg_vel(s);
487 0 : double velX = get_avg_velX(s);
488 0 : file << _counter << " ; " << vel << " ; " << velX << std::endl;
489 0 : file.close();
490 0 : }
491 :
492 : /// ------------------------------------------------------------------------------------------------------
493 : /// End of Pint methods ------ End of Pint methods ------ End of Pint methods ------ End of Pint methods
494 : /// ------------------------------------------------------------------------------------------------------
495 :
496 : /** calls stream() and collide() and swaps the fields
497 : * @brief collide-stream algorithm for the Lattice-Boltzmann method */
498 51 : void collidestream() {
499 51 : #pragma omp parallel for
500 : for (int z = 1; z < _domainSizeZ + 1; z++) {
501 : for (int y = 1; y < _domainSizeY + 1; y++) {
502 : for (int x = 1; x < _domainSizeX + 1; x++) {
503 : const int index = get(x, y, z);
504 : if (_flag[index] == FLUID) {
505 : stream(index);
506 : collide(index, x, y, z);
507 : }
508 : }
509 : }
510 : }
511 : // swap fields
512 51 : double* swap = _pdf1;
513 51 : _pdf1 = _pdf2;
514 51 : _pdf2 = swap;
515 51 : }
516 :
517 : /** @brief the stream part of the LB algorithm (from pdf1 to pdf2) */
518 408000 : void stream(int index) {
519 408000 : const int pI = 19 * index;
520 4080000 : for (int q = 0; q < 9; q++) {
521 3672000 : const int nb = 19 * (_C[q][0] + _C[q][1] * _xO + _C[q][2] * _yO);
522 3672000 : _pdf2[pI + q] = _pdf1[pI + q - nb];
523 3672000 : _pdf2[pI + 18 - q] = _pdf1[pI + 18 - q + nb];
524 : }
525 408000 : _pdf2[pI + 9] = _pdf1[pI + 9];
526 408000 : }
527 :
528 : /** @brieff the collide step within pdf2 */
529 408000 : void collide(int index, int x, int y, int z) {
530 : // index of start of cell-local pdfs in AoS
531 408000 : const int pI = 19 * index;
532 : // compute and store density, velocity
533 408000 : double* vel = &_vel[3 * index];
534 408000 : computeDensityAndVelocity(vel, _density[index], &_pdf2[pI]);
535 : // collide (BGK); always handle pdfs no. q and inv(q)=18-q in one step
536 408000 : const double u2 = 1.0 - 1.5 * (vel[0] * vel[0] + vel[1] * vel[1] + vel[2] * vel[2]);
537 : // pdf 0,18
538 408000 : double cu = -vel[1] - vel[2];
539 408000 : int nb = -_xO - _yO;
540 408000 : double feq = _W[0] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
541 408000 : _pdf2[pI] -= _omega * (_pdf2[pI] - feq);
542 408000 : boundary(_pdf2, pI, x, y, z, 0, _flag[index + nb], pI + 19 * nb);
543 408000 : feq -= 6.0 * _W[0] * _density[index] * cu;
544 408000 : _pdf2[pI + 18] -= _omega * (_pdf2[pI + 18] - feq);
545 408000 : boundary(_pdf2, pI, x, y, z, 18, _flag[index - nb], pI - 19 * nb);
546 : // pdf 1,17
547 408000 : cu = -vel[0] - vel[2];
548 408000 : nb = -1 - _yO;
549 408000 : feq = _W[1] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
550 408000 : _pdf2[pI + 1] -= _omega * (_pdf2[pI + 1] - feq);
551 408000 : boundary(_pdf2, pI, x, y, z, 1, _flag[index + nb], pI + 19 * nb);
552 408000 : feq -= 6.0 * _W[1] * _density[index] * cu;
553 408000 : _pdf2[pI + 17] -= _omega * (_pdf2[pI + 17] - feq);
554 408000 : boundary(_pdf2, pI, x, y, z, 17, _flag[index - nb], pI - 19 * nb);
555 : // pdf 2,16
556 408000 : cu = -vel[2];
557 408000 : nb = -_yO;
558 408000 : feq = _W[2] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
559 408000 : _pdf2[pI + 2] -= _omega * (_pdf2[pI + 2] - feq);
560 408000 : boundary(_pdf2, pI, x, y, z, 2, _flag[index + nb], pI + 19 * nb);
561 408000 : feq -= 6.0 * _W[2] * _density[index] * cu;
562 408000 : _pdf2[pI + 16] -= _omega * (_pdf2[pI + 16] - feq);
563 408000 : boundary(_pdf2, pI, x, y, z, 16, _flag[index - nb], pI - 19 * nb);
564 : // pdf 3,15
565 408000 : cu = vel[0] - vel[2];
566 408000 : nb = 1 - _yO;
567 408000 : feq = _W[3] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
568 408000 : _pdf2[pI + 3] -= _omega * (_pdf2[pI + 3] - feq);
569 408000 : boundary(_pdf2, pI, x, y, z, 3, _flag[index + nb], pI + 19 * nb);
570 408000 : feq -= 6.0 * _W[3] * _density[index] * cu;
571 408000 : _pdf2[pI + 15] -= _omega * (_pdf2[pI + 15] - feq);
572 408000 : boundary(_pdf2, pI, x, y, z, 15, _flag[index - nb], pI - 19 * nb);
573 : // pdf 4,14
574 408000 : cu = vel[1] - vel[2];
575 408000 : nb = _xO - _yO;
576 408000 : feq = _W[4] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
577 408000 : _pdf2[pI + 4] -= _omega * (_pdf2[pI + 4] - feq);
578 408000 : boundary(_pdf2, pI, x, y, z, 4, _flag[index + nb], pI + 19 * nb);
579 408000 : feq -= 6.0 * _W[4] * _density[index] * cu;
580 408000 : _pdf2[pI + 14] -= _omega * (_pdf2[pI + 14] - feq);
581 408000 : boundary(_pdf2, pI, x, y, z, 14, _flag[index - nb], pI - 19 * nb);
582 : // pdf 5,13
583 408000 : cu = -vel[0] - vel[1];
584 408000 : nb = -1 - _xO;
585 408000 : feq = _W[5] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
586 408000 : _pdf2[pI + 5] -= _omega * (_pdf2[pI + 5] - feq);
587 408000 : boundary(_pdf2, pI, x, y, z, 5, _flag[index + nb], pI + 19 * nb);
588 408000 : feq -= 6.0 * _W[5] * _density[index] * cu;
589 408000 : _pdf2[pI + 13] -= _omega * (_pdf2[pI + 13] - feq);
590 408000 : boundary(_pdf2, pI, x, y, z, 13, _flag[index - nb], pI - 19 * nb);
591 : // pdf 6,12
592 408000 : cu = -vel[1];
593 408000 : nb = -_xO;
594 408000 : feq = _W[6] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
595 408000 : _pdf2[pI + 6] -= _omega * (_pdf2[pI + 6] - feq);
596 408000 : boundary(_pdf2, pI, x, y, z, 6, _flag[index + nb], pI + 19 * nb);
597 408000 : feq -= 6.0 * _W[6] * _density[index] * cu;
598 408000 : _pdf2[pI + 12] -= _omega * (_pdf2[pI + 12] - feq);
599 408000 : boundary(_pdf2, pI, x, y, z, 12, _flag[index - nb], pI - 19 * nb);
600 : // pdf 7,11
601 408000 : cu = vel[0] - vel[1];
602 408000 : nb = 1 - _xO;
603 408000 : feq = _W[7] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
604 408000 : _pdf2[pI + 7] -= _omega * (_pdf2[pI + 7] - feq);
605 408000 : boundary(_pdf2, pI, x, y, z, 7, _flag[index + nb], pI + 19 * nb);
606 408000 : feq -= 6.0 * _W[7] * _density[index] * cu;
607 408000 : _pdf2[pI + 11] -= _omega * (_pdf2[pI + 11] - feq);
608 408000 : boundary(_pdf2, pI, x, y, z, 11, _flag[index - nb], pI - 19 * nb);
609 : // pdf 8,10
610 408000 : cu = -vel[0];
611 408000 : nb = -1;
612 408000 : feq = _W[8] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
613 408000 : _pdf2[pI + 8] -= _omega * (_pdf2[pI + 8] - feq);
614 408000 : boundary(_pdf2, pI, x, y, z, 8, _flag[index + nb], pI + 19 * nb);
615 408000 : feq -= 6.0 * _W[8] * _density[index] * cu;
616 408000 : _pdf2[pI + 10] -= _omega * (_pdf2[pI + 10] - feq);
617 408000 : boundary(_pdf2, pI, x, y, z, 10, _flag[index - nb], pI - 19 * nb);
618 : // pdf 9
619 408000 : _pdf2[pI + 9] -= _omega * (_pdf2[pI + 9] - _W[9] * _density[index] * u2);
620 408000 : }
621 :
622 : /** @brief takes care of the correct boundary treatment for the LB method
623 : * @param pdf particle distribution function
624 : * @param index start index for current cell in pdf-array
625 : * @param x the position in x direction of the cell
626 : * @param y the position in y direction of the cell
627 : * @param z the position in z direction of the cell
628 : * @param q distribution function number
629 : * @param flag boundary flag of neighbouring cell
630 : * @param nbIndex index of neighbouring cell */
631 7344000 : void boundary(double* const pdf, int index, int x, int y, int z, int q, const Flag& flag, int nbIndex) {
632 7344000 : if (flag != FLUID) {
633 599760 : if (flag == NO_SLIP) {
634 : // half-way bounce back
635 102000 : pdf[nbIndex + 18 - q] = pdf[index + q];
636 497760 : } else if (flag == MOVING_WALL) {
637 : // half-way bounce back + moving wall acceleration (only x-direction for
638 : // wall supported at the moment)
639 102000 : pdf[nbIndex + 18 - q] =
640 102000 : pdf[index + q] - 6.0 * _W[q] * _density[index / 19] * (_C[q][0] * _wallVelocity[0] + _C[q][1] * _wallVelocity[1] + _C[q][2] * _wallVelocity[2]);
641 395760 : } else if (flag == PERIODIC) {
642 : // periodic treatment
643 0 : int target[3] = {x, y, z};
644 0 : if (target[0] + _C[q][0] == 0) {
645 0 : target[0] = _domainSizeX + 1;
646 0 : } else if (target[0] + _C[q][0] == _domainSizeX + 1) {
647 0 : target[0] = 0;
648 : }
649 0 : if (target[1] + _C[q][1] == 0) {
650 0 : target[1] = _domainSizeY + 1;
651 0 : } else if (target[1] + _C[q][1] == _domainSizeY + 1) {
652 0 : target[1] = 0;
653 : }
654 0 : if (target[2] + _C[q][2] == 0) {
655 0 : target[2] = _domainSizeZ + 1;
656 0 : } else if (target[2] + _C[q][2] == _domainSizeZ + 1) {
657 0 : target[2] = 0;
658 : }
659 0 : const int periodicNb = target[0] + (_domainSizeX + 2) * (target[1] + (_domainSizeY + 2) * target[2]);
660 0 : pdf[19 * periodicNb + q] = pdf[index + q];
661 : }
662 : }
663 7344000 : }
664 :
665 : /** @brief refers to the LB method; computes density and velocity on pdf
666 : * @param vel velocity
667 : * @param density density
668 : * @param pdf partial distribution function */
669 965240 : static void computeDensityAndVelocity(double* const vel, double& density, const double* const pdf) {
670 965240 : vel[0] = -(pdf[1] + pdf[5] + pdf[8] + pdf[11] + pdf[15]);
671 965240 : density = pdf[3] + pdf[7] + pdf[10] + pdf[13] + pdf[17];
672 965240 : vel[1] = (pdf[4] + pdf[11] + pdf[12] + pdf[13] + pdf[18]) - (pdf[0] + pdf[5] + pdf[6] + pdf[7] + pdf[14]);
673 965240 : vel[0] = density + vel[0];
674 965240 : density = density + pdf[0] + pdf[1] + pdf[2] + pdf[4] + pdf[5] + pdf[6] + pdf[8] + pdf[9] + pdf[11] + pdf[12] + pdf[14] + pdf[15] + pdf[16] + pdf[18];
675 965240 : vel[2] = (pdf[14] + pdf[15] + pdf[16] + pdf[17] + pdf[18]) - (pdf[0] + pdf[1] + pdf[2] + pdf[3] + pdf[4]);
676 965240 : vel[0] = vel[0] / density;
677 965240 : vel[1] = vel[1] / density;
678 965240 : vel[2] = vel[2] / density;
679 965240 : }
680 :
681 : /** takes care of communication across one face in one direction.
682 : * @param pdf partial distribution function
683 : * @param sendBuffer send buffer
684 : * @param recvBuffer receive buffer
685 : * @param nbFlagTo direction into which message is sent
686 : * @param nbFlagFrom direction from which message is received
687 : * @param startSend 3d coordinates that define the start of the data to be
688 : * sent to neighbouring process
689 : * @param endSend 3d coordinates that define the end of the data to to be
690 : * sent to neighbouring process
691 : * @param startRecv 3d coordinates that define the start of the data to be
692 : * received from neighbouring process
693 : * @param endRecv 3d coordinates that define the end of the data to be
694 : * received from neighbouring process */
695 306 : void communicatePart(double* pdf, double* sendBuffer, double* recvBuffer, NbFlag nbFlagTo, NbFlag nbFlagFrom, tarch::la::Vector<3, int> startSend,
696 : tarch::la::Vector<3, int> endSend, tarch::la::Vector<3, int> startRecv, tarch::la::Vector<3, int> endRecv) {
697 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
698 : // directions that point to LEFT/RIGHT,... -> same ordering as enums!
699 306 : const int directions[6][5] = {{1, 5, 8, 11, 15}, {3, 7, 10, 13, 17}, {4, 11, 12, 13, 18}, {0, 5, 6, 7, 14}, {0, 1, 2, 3, 4}, {14, 15, 16, 17, 18}};
700 306 : MPI_Request requests[2];
701 306 : MPI_Status status[2];
702 306 : tarch::la::Vector<2, int> plane;
703 306 : tarch::la::Vector<2, int> domainSize;
704 : // find out plane coordinates
705 306 : if (nbFlagTo == LEFT || nbFlagTo == RIGHT) {
706 102 : plane[0] = 1;
707 102 : plane[1] = 2;
708 102 : domainSize[0] = _domainSizeY;
709 102 : domainSize[1] = _domainSizeZ;
710 204 : } else if (nbFlagTo == FRONT || nbFlagTo == BACK) {
711 102 : plane[0] = 0;
712 102 : plane[1] = 2;
713 102 : domainSize[0] = _domainSizeX;
714 102 : domainSize[1] = _domainSizeZ;
715 102 : } else if (nbFlagTo == TOP || nbFlagTo == BOTTOM) {
716 102 : plane[0] = 0;
717 102 : plane[1] = 1;
718 102 : domainSize[0] = _domainSizeX;
719 102 : domainSize[1] = _domainSizeY;
720 : } else {
721 0 : std::cout << "ERROR LBCouetteSolver::communicatePart: d >2 or d < 0!" << std::endl;
722 0 : exit(EXIT_FAILURE);
723 : }
724 : // extract data and write to send buffer
725 306 : tarch::la::Vector<3, int> coords(0);
726 4488 : for (coords[2] = startSend[2]; coords[2] < endSend[2]; coords[2]++) {
727 49266 : for (coords[1] = startSend[1]; coords[1] < endSend[1]; coords[1]++) {
728 180132 : for (coords[0] = startSend[0]; coords[0] < endSend[0]; coords[0]++) {
729 810288 : for (int q = 0; q < 5; q++) {
730 675240 : sendBuffer[q + 5 * getParBuf(coords[plane[0]], coords[plane[1]], domainSize[0], domainSize[1])] =
731 675240 : pdf[directions[nbFlagTo][q] + 19 * get(coords[0], coords[1], coords[2])];
732 : }
733 : }
734 : }
735 : }
736 : // send and receive data
737 306 : MPI_Irecv(recvBuffer, (domainSize[0] + 2) * (domainSize[1] + 2) * 5, MPI_DOUBLE, _parallelNeighbours[nbFlagFrom], 1000,
738 306 : coupling::indexing::IndexingService<3>::getInstance().getComm(), &requests[0]);
739 306 : MPI_Isend(sendBuffer, (domainSize[0] + 2) * (domainSize[1] + 2) * 5, MPI_DOUBLE, _parallelNeighbours[nbFlagTo], 1000,
740 306 : coupling::indexing::IndexingService<3>::getInstance().getComm(), &requests[1]);
741 306 : MPI_Waitall(2, requests, status);
742 : // write data back to pdf field
743 306 : if (_parallelNeighbours[nbFlagFrom] != MPI_PROC_NULL) {
744 4284 : for (coords[2] = startRecv[2]; coords[2] < endRecv[2]; coords[2]++) {
745 46920 : for (coords[1] = startRecv[1]; coords[1] < endRecv[1]; coords[1]++) {
746 128520 : for (coords[0] = startRecv[0]; coords[0] < endRecv[0]; coords[0]++) {
747 514080 : for (int q = 0; q < 5; q++) {
748 428400 : if (_flag[get(coords[0], coords[1], coords[2])] == PARALLEL_BOUNDARY) {
749 428400 : pdf[directions[nbFlagTo][q] + 19 * get(coords[0], coords[1], coords[2])] =
750 428400 : recvBuffer[q + 5 * getParBuf(coords[plane[0]], coords[plane[1]], domainSize[0], domainSize[1])];
751 : }
752 : }
753 : }
754 : }
755 : }
756 : }
757 : #endif
758 306 : }
759 :
760 : /** @brief comunicates the boundary field data between the different processes
761 : */
762 51 : void communicate() {
763 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
764 : // send from right to left
765 102 : communicatePart(_pdf1, _sendBufferX, _recvBufferX, LEFT, RIGHT, tarch::la::Vector<3, int>(1, 1, 1),
766 51 : tarch::la::Vector<3, int>(2, _domainSizeY + 1, _domainSizeZ + 1), tarch::la::Vector<3, int>(_domainSizeX + 1, 1, 1),
767 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 1, _domainSizeZ + 1));
768 : // send from left to right
769 102 : communicatePart(_pdf1, _sendBufferX, _recvBufferX, RIGHT, LEFT, tarch::la::Vector<3, int>(_domainSizeX, 1, 1),
770 51 : tarch::la::Vector<3, int>(_domainSizeX + 1, _domainSizeY + 1, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, 1, 1),
771 51 : tarch::la::Vector<3, int>(1, _domainSizeY + 1, _domainSizeZ + 1));
772 : // send from back to front
773 102 : communicatePart(_pdf1, _sendBufferY, _recvBufferY, FRONT, BACK, tarch::la::Vector<3, int>(0, 1, 1),
774 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, 2, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, _domainSizeY + 1, 1),
775 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, _domainSizeZ + 1));
776 : // send from front to back
777 102 : communicatePart(_pdf1, _sendBufferY, _recvBufferY, BACK, FRONT, tarch::la::Vector<3, int>(0, _domainSizeY, 1),
778 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 1, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, 0, 1),
779 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, 1, _domainSizeZ + 1));
780 : // send from top to bottom
781 102 : communicatePart(_pdf1, _sendBufferZ, _recvBufferZ, BOTTOM, TOP, tarch::la::Vector<3, int>(0, 0, 1),
782 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, 2), tarch::la::Vector<3, int>(0, 0, _domainSizeZ + 1),
783 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, _domainSizeZ + 2));
784 : // send from bottom to top
785 102 : communicatePart(_pdf1, _sendBufferZ, _recvBufferZ, TOP, BOTTOM, tarch::la::Vector<3, int>(0, 0, _domainSizeZ),
786 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, 0, 0),
787 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, 1));
788 : #endif
789 51 : }
790 :
791 : /** @brief relaxation frequency */
792 : const double _omega;
793 : /** @brief velocity of moving wall of Couette flow */
794 : tarch::la::Vector<3, double> _wallVelocity;
795 : int _pdfsize{0};
796 : /** @brief partical distribution function field */
797 : double* _pdf1{NULL};
798 : /** @brief partial distribution function field (stores the old time step)*/
799 : double* _pdf2{NULL};
800 : /** @brief lattice velocities*/
801 : const int _C[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},
802 : {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}};
803 : /** @brief lattice weights */
804 : const double _W[19]{1.0 / 36.0, 1.0 / 36.0, 1.0 / 18.0, 1.0 / 36.0, 1.0 / 36.0, 1.0 / 36.0, 1.0 / 18.0, 1.0 / 36.0, 1.0 / 18.0, 1.0 / 3.0,
805 : 1.0 / 18.0, 1.0 / 36.0, 1.0 / 18.0, 1.0 / 36.0, 1.0 / 36.0, 1.0 / 36.0, 1.0 / 18.0, 1.0 / 36.0, 1.0 / 36.0};
806 : /** @brief enables avg_vel CSV output */
807 : const bool _plotAverageVelocity;
808 : };
809 :
810 : #endif // _MOLECULARDYNAMICS_COUPLING_SOLVERS_LBCOUETTESOLVER_H_
|
