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