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 4860 : LBCouetteSolverState(int size) : _pdf(size, 0) {}
24 :
25 152 : 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 46 : 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 348 : _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 : 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 : #pragma omp parallel for
114 308821 : for (int i = 0; i < (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2); i++) {
115 6175840 : for (int q = 0; q < 19; q++) {
116 5867048 : _pdf1[get(i) * 19 + q] = _W[q];
117 5867048 : _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 :
210 : // loop over all received cells
211 0 : for (auto pair : md2macroBuffer) {
212 0 : I01 idx;
213 0 : const coupling::datastructures::CouplingCell<3>* couplingCell;
214 0 : std::tie(couplingCell, idx) = pair;
215 : // determine cell index of this cell in LB domain
216 0 : tarch::la::Vector<3, unsigned int> globalCellCoords{idx.get()};
217 0 : globalCellCoords[0] = (globalCellCoords[0] + _offset[0]) - _coords[0] * _avgDomainSizeX;
218 0 : globalCellCoords[1] = (globalCellCoords[1] + _offset[1]) - _coords[1] * _avgDomainSizeY;
219 0 : globalCellCoords[2] = (globalCellCoords[2] + _offset[2]) - _coords[2] * _avgDomainSizeZ;
220 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
221 : std::cout << "Process coords: " << _coords << ": GlobalCellCoords for index " << idx << ": " << globalCellCoords << std::endl;
222 : #endif
223 0 : const int index = get(globalCellCoords[0], globalCellCoords[1], globalCellCoords[2]);
224 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
225 0 : if (_flag[index] != MD_BOUNDARY) {
226 0 : std::cout << "ERROR LBCouetteSolver::setMDBoundaryValues(): Cell " << index << " is no MD boundary cell!" << std::endl;
227 0 : exit(EXIT_FAILURE);
228 : }
229 : #endif
230 : // set velocity value and pdfs in MD boundary cell (before streaming); the
231 : // boundary velocities are interpolated between the neighbouring and this
232 : // cell. This interpolation is valid for FLUID-MD_BOUNDARY neighbouring
233 : // relations only. determine local velocity received from MaMiCo and
234 : // convert it to LB units; store the velocity in _vel
235 0 : tarch::la::Vector<3, double> localVel((1.0 / couplingCell->getMacroscopicMass()) * (_dt / _dx) * couplingCell->getMacroscopicMomentum());
236 0 : for (unsigned int d = 0; d < 3; d++) {
237 0 : _vel[3 * index + d] = localVel[d];
238 : }
239 : // loop over all pdfs and set them according to interpolated moving-wall
240 : // conditions
241 0 : for (unsigned int q = 0; q < 19; q++) {
242 : // index of neighbour cell; only if cell is located inside local domain
243 0 : if (((int)globalCellCoords[0] + _C[q][0] > 0) && ((int)globalCellCoords[0] + _C[q][0] < _domainSizeX + 1) &&
244 0 : ((int)globalCellCoords[1] + _C[q][1] > 0) && ((int)globalCellCoords[1] + _C[q][1] < _domainSizeY + 1) &&
245 0 : ((int)globalCellCoords[2] + _C[q][2] > 0) && ((int)globalCellCoords[2] + _C[q][2] < _domainSizeZ + 1)) {
246 0 : const int nbIndex = get((_C[q][0] + globalCellCoords[0]), (_C[q][1] + globalCellCoords[1]), (_C[q][2] + globalCellCoords[2]));
247 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]),
248 0 : 0.5 * (_vel[3 * index + 2] + _vel[3 * nbIndex + 2]));
249 0 : _pdf1[19 * index + q] =
250 0 : _pdf1[19 * nbIndex + 18 - q] -
251 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]);
252 : }
253 : }
254 : }
255 : }
256 :
257 : /** @brief returns velocity at a certain position
258 : * @param pos position for which the velocity will be returned
259 : * @returns the velocity vector for the position */
260 300 : tarch::la::Vector<3, double> getVelocity(tarch::la::Vector<3, double> pos) const override {
261 300 : tarch::la::Vector<3, unsigned int> coords;
262 300 : const tarch::la::Vector<3, double> domainOffset(_coords[0] * _dx * _avgDomainSizeX, _coords[1] * _dx * _avgDomainSizeY, _coords[2] * _dx * _avgDomainSizeZ);
263 : // check pos-data for process locality (todo: put this in debug mode in
264 : // future releases)
265 300 : if ((pos[0] < domainOffset[0]) || (pos[0] > domainOffset[0] + _domainSizeX * _dx) || (pos[1] < domainOffset[1]) ||
266 600 : (pos[1] > domainOffset[1] + _domainSizeY * _dx) || (pos[2] < domainOffset[2]) || (pos[2] > domainOffset[2] + _domainSizeZ * _dx)) {
267 0 : std::cout << "ERROR LBCouetteSolver::getVelocity(): Position " << pos << " out of range!" << std::endl;
268 0 : std::cout << "domainOffset = " << domainOffset << std::endl;
269 0 : std::cout << "_domainSizeX = " << _domainSizeX << std::endl;
270 0 : std::cout << "_domainSizeY = " << _domainSizeY << std::endl;
271 0 : std::cout << "_domainSizeZ = " << _domainSizeZ << std::endl;
272 0 : std::cout << "_dx = " << _dx << std::endl;
273 0 : exit(EXIT_FAILURE);
274 : }
275 : // compute index for respective cell (_dx+... for ghost cells); use coords
276 : // to store local cell coordinates
277 1200 : for (unsigned int d = 0; d < 3; d++) {
278 900 : coords[d] = (unsigned int)((_dx + pos[d] - domainOffset[d]) / _dx);
279 : }
280 300 : const int index = get(coords[0], coords[1], coords[2]);
281 300 : tarch::la::Vector<3, double> vel(0.0);
282 : // extract and scale velocity to "real"=MD units
283 1200 : for (int d = 0; d < 3; d++) {
284 900 : vel[d] = _dx / _dt * _vel[3 * index + d];
285 : }
286 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
287 : std::cout << "Position " << pos << " corresponds to cell: " << coords << "; vel=" << vel << std::endl;
288 : #endif
289 300 : return vel;
290 : }
291 :
292 : /** @brief returns density at a certain position
293 : * @param pos position for which the density will be returned
294 : * @returns the density vector for the position */
295 300 : double getDensity(tarch::la::Vector<3, double> pos) const override {
296 300 : tarch::la::Vector<3, unsigned int> coords;
297 300 : const tarch::la::Vector<3, double> domainOffset(_coords[0] * _dx * _avgDomainSizeX, _coords[1] * _dx * _avgDomainSizeY, _coords[2] * _dx * _avgDomainSizeZ);
298 : // check pos-data for process locality (todo: put this in debug mode in
299 : // future releases)
300 300 : if ((pos[0] < domainOffset[0]) || (pos[0] > domainOffset[0] + _domainSizeX * _dx) || (pos[1] < domainOffset[1]) ||
301 600 : (pos[1] > domainOffset[1] + _domainSizeY * _dx) || (pos[2] < domainOffset[2]) || (pos[2] > domainOffset[2] + _domainSizeZ * _dx)) {
302 0 : std::cout << "ERROR LBCouetteSolver::getDensity(): Position " << pos << " out of range!" << std::endl;
303 0 : exit(EXIT_FAILURE);
304 : }
305 : // compute index for respective cell (_dx+... for ghost cells); use coords
306 : // to store local cell coordinates
307 1200 : for (unsigned int d = 0; d < 3; d++) {
308 900 : coords[d] = (unsigned int)((_dx + pos[d] - domainOffset[d]) / _dx);
309 : }
310 300 : const int index = get(coords[0], coords[1], coords[2]);
311 300 : return _density[index];
312 : }
313 :
314 : /** @brief changes the velocity at the moving wall (z=0)
315 : * @param wallVelocity the velocity will be set at the moving wall */
316 0 : virtual void setWallVelocity(const tarch::la::Vector<3, double> wallVelocity) override { _wallVelocity = (_dt / _dx) * wallVelocity; }
317 :
318 : /// ------------------------------------------------------------------------------------------------------
319 : /// Pint methods ------ Pint methods ------ Pint methods ------ Pint methods ------ Pint methods
320 : /// ------------------------------------------------------------------------------------------------------
321 :
322 80 : std::unique_ptr<State> getState() override {
323 80 : computeDensityAndVelocityEverywhere();
324 80 : return std::make_unique<LBCouetteSolverState>(_pdfsize, _pdf1);
325 : }
326 :
327 23 : void setState(const std::unique_ptr<State>& input, int cycle) override {
328 23 : const LBCouetteSolverState* state = dynamic_cast<const LBCouetteSolverState*>(input.get());
329 :
330 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
331 23 : if (state == nullptr) {
332 0 : std::cout << "ERROR LBCouetteSolver setState() wrong state type" << std::endl;
333 0 : exit(EXIT_FAILURE);
334 : }
335 : #endif
336 :
337 23 : std::copy(state->getData(), state->getData() + _pdfsize, _pdf1);
338 23 : computeDensityAndVelocityEverywhere();
339 :
340 23 : _counter = cycle;
341 23 : }
342 :
343 14 : std::unique_ptr<State> operator()(const std::unique_ptr<State>& input, int cycle) override {
344 14 : setState(input, cycle);
345 :
346 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
347 14 : if (_mode != Mode::supervising) {
348 0 : std::cout << "ERROR LBCouetteSolver operator() called but not in supervising mode" << std::endl;
349 0 : exit(EXIT_FAILURE);
350 : }
351 : #endif
352 :
353 14 : advance(_dt_pint);
354 14 : return getState();
355 : }
356 :
357 12 : Mode getMode() const override { return _mode; }
358 :
359 : /**
360 : * This will create a new instance of this LBCouetteSolver
361 : * In the supervisor, setMDBoundary() has not been called, independently of state of couette solver in coupling mode
362 : * The supervisor can run with a modified viscosity
363 : * The supervisor's operator() advances many coupling cycles at once, interval is stored in _dt_pint
364 : */
365 72 : std::unique_ptr<PintableMacroSolver> getSupervisor(int num_cycles, double visc_multiplier) const override {
366 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
367 72 : if (_mode == Mode::supervising) {
368 0 : std::cout << "ERROR LBCouetteSolver getSupervisor(): already in supervising mode" << std::endl;
369 0 : exit(EXIT_FAILURE);
370 : }
371 : #endif
372 :
373 72 : int numThreads = 1;
374 : #if defined(_OPENMP)
375 : numThreads = omp_get_num_threads();
376 : #endif
377 :
378 144 : auto res = std::make_unique<LBCouetteSolver>(_channelheight, _wallVelocity * _dx / _dt, _kinVisc * visc_multiplier, _dx, _dt, _plotEveryTimestep, _plotAverageVelocity,
379 216 : _filestem + std::string("_supervising"), _processes, numThreads, _scen);
380 :
381 72 : res->_mode = Mode::supervising;
382 72 : res->_dt_pint = _dt * num_cycles;
383 :
384 144 : return res;
385 72 : }
386 :
387 0 : void print(std::ostream& os) const override {
388 0 : if (_mode == Mode::supervising)
389 0 : os << "<LBCouetteSolver instance in supervising mode >";
390 0 : if (_mode == Mode::coupling)
391 0 : os << "<LBCouetteSolver instance in coupling mode >";
392 0 : }
393 :
394 29 : double get_avg_vel(const std::unique_ptr<State>& state) const override {
395 29 : double vel[3];
396 29 : double density;
397 29 : double res[3]{0, 0, 0};
398 53269 : for (int i = 0; i < _pdfsize; i += 19) {
399 53240 : LBCouetteSolver::computeDensityAndVelocity(vel, density, state->getData() + i);
400 53240 : res[0] += vel[0];
401 53240 : res[1] += vel[1];
402 53240 : res[2] += vel[2];
403 : }
404 29 : if (_pdfsize > 0) {
405 5 : res[0] /= (_pdfsize / 19);
406 5 : res[1] /= (_pdfsize / 19);
407 5 : res[2] /= (_pdfsize / 19);
408 : }
409 29 : return std::sqrt(res[0] * res[0] + res[1] * res[1] + res[2] * res[2]);
410 : }
411 :
412 : private:
413 : Mode _mode;
414 : double _dt_pint;
415 :
416 132 : void computeDensityAndVelocityEverywhere() {
417 132 : if (skipRank())
418 : return;
419 1323 : for (int z = 1; z < _domainSizeZ + 1; z++) {
420 26460 : for (int y = 1; y < _domainSizeY + 1; y++) {
421 529200 : for (int x = 1; x < _domainSizeX + 1; x++) {
422 504000 : const int index = get(x, y, z);
423 504000 : const int pI = 19 * index;
424 504000 : double* vel = &_vel[3 * index];
425 504000 : computeDensityAndVelocity(vel, _density[index], &_pdf1[pI]);
426 : }
427 : }
428 : }
429 : }
430 :
431 51 : void plot_avg_vel(){
432 51 : if (!_plotAverageVelocity) return;
433 :
434 0 : int rank = 0;
435 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
436 0 : MPI_Comm_rank(coupling::indexing::IndexingService<3>::getInstance().getComm(), &rank);
437 : #endif
438 0 : std::stringstream ss;
439 0 : ss << _filestem << "_r" << rank;
440 0 : if (_scen != nullptr) {
441 0 : auto ts = _scen->getTimeIntegrationService();
442 0 : if (ts != nullptr) {
443 0 : if (ts->isPintEnabled())
444 0 : ss << "_i" << ts->getInteration();
445 : }
446 : }
447 0 : ss << ".csv";
448 0 : std::string filename = ss.str();
449 0 : std::ofstream file(filename.c_str(), _counter==0 ? std::ofstream::out : std::ofstream::app);
450 0 : if (!file.is_open()) {
451 0 : std::cout << "ERROR LBCouetteSolver::plot_avg_vel(): Could not open file " << filename << "!" << std::endl;
452 0 : exit(EXIT_FAILURE);
453 : }
454 :
455 0 : if(_counter == 0){
456 0 : file << "coupling_cycle ; avg_vel" << std::endl;
457 : }
458 :
459 0 : std::unique_ptr<State> s = std::make_unique<LBCouetteSolverState>(_pdfsize, _pdf1);
460 0 : double vel = get_avg_vel(s);
461 0 : file << _counter << " ; " << vel << std::endl;
462 0 : file.close();
463 0 : }
464 :
465 : /// ------------------------------------------------------------------------------------------------------
466 : /// End of Pint methods ------ End of Pint methods ------ End of Pint methods ------ End of Pint methods
467 : /// ------------------------------------------------------------------------------------------------------
468 :
469 : /** calls stream() and collide() and swaps the fields
470 : * @brief collide-stream algorithm for the Lattice-Boltzmann method */
471 51 : void collidestream() {
472 : #pragma omp parallel for
473 1071 : for (int z = 1; z < _domainSizeZ + 1; z++) {
474 21420 : for (int y = 1; y < _domainSizeY + 1; y++) {
475 428400 : for (int x = 1; x < _domainSizeX + 1; x++) {
476 408000 : const int index = get(x, y, z);
477 408000 : if (_flag[index] == FLUID) {
478 408000 : stream(index);
479 408000 : collide(index, x, y, z);
480 : }
481 : }
482 : }
483 : }
484 : // swap fields
485 51 : double* swap = _pdf1;
486 51 : _pdf1 = _pdf2;
487 51 : _pdf2 = swap;
488 51 : }
489 :
490 : /** @brief the stream part of the LB algorithm (from pdf1 to pdf2) */
491 408000 : void stream(int index) {
492 408000 : const int pI = 19 * index;
493 4080000 : for (int q = 0; q < 9; q++) {
494 3672000 : const int nb = 19 * (_C[q][0] + _C[q][1] * _xO + _C[q][2] * _yO);
495 3672000 : _pdf2[pI + q] = _pdf1[pI + q - nb];
496 3672000 : _pdf2[pI + 18 - q] = _pdf1[pI + 18 - q + nb];
497 : }
498 408000 : _pdf2[pI + 9] = _pdf1[pI + 9];
499 408000 : }
500 :
501 : /** @brieff the collide step within pdf2 */
502 408000 : void collide(int index, int x, int y, int z) {
503 : // index of start of cell-local pdfs in AoS
504 408000 : const int pI = 19 * index;
505 : // compute and store density, velocity
506 408000 : double* vel = &_vel[3 * index];
507 408000 : computeDensityAndVelocity(vel, _density[index], &_pdf2[pI]);
508 : // collide (BGK); always handle pdfs no. q and inv(q)=18-q in one step
509 408000 : const double u2 = 1.0 - 1.5 * (vel[0] * vel[0] + vel[1] * vel[1] + vel[2] * vel[2]);
510 : // pdf 0,18
511 408000 : double cu = -vel[1] - vel[2];
512 408000 : int nb = -_xO - _yO;
513 408000 : double feq = _W[0] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
514 408000 : _pdf2[pI] -= _omega * (_pdf2[pI] - feq);
515 408000 : boundary(_pdf2, pI, x, y, z, 0, _flag[index + nb], pI + 19 * nb);
516 408000 : feq -= 6.0 * _W[0] * _density[index] * cu;
517 408000 : _pdf2[pI + 18] -= _omega * (_pdf2[pI + 18] - feq);
518 408000 : boundary(_pdf2, pI, x, y, z, 18, _flag[index - nb], pI - 19 * nb);
519 : // pdf 1,17
520 408000 : cu = -vel[0] - vel[2];
521 408000 : nb = -1 - _yO;
522 408000 : feq = _W[1] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
523 408000 : _pdf2[pI + 1] -= _omega * (_pdf2[pI + 1] - feq);
524 408000 : boundary(_pdf2, pI, x, y, z, 1, _flag[index + nb], pI + 19 * nb);
525 408000 : feq -= 6.0 * _W[1] * _density[index] * cu;
526 408000 : _pdf2[pI + 17] -= _omega * (_pdf2[pI + 17] - feq);
527 408000 : boundary(_pdf2, pI, x, y, z, 17, _flag[index - nb], pI - 19 * nb);
528 : // pdf 2,16
529 408000 : cu = -vel[2];
530 408000 : nb = -_yO;
531 408000 : feq = _W[2] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
532 408000 : _pdf2[pI + 2] -= _omega * (_pdf2[pI + 2] - feq);
533 408000 : boundary(_pdf2, pI, x, y, z, 2, _flag[index + nb], pI + 19 * nb);
534 408000 : feq -= 6.0 * _W[2] * _density[index] * cu;
535 408000 : _pdf2[pI + 16] -= _omega * (_pdf2[pI + 16] - feq);
536 408000 : boundary(_pdf2, pI, x, y, z, 16, _flag[index - nb], pI - 19 * nb);
537 : // pdf 3,15
538 408000 : cu = vel[0] - vel[2];
539 408000 : nb = 1 - _yO;
540 408000 : feq = _W[3] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
541 408000 : _pdf2[pI + 3] -= _omega * (_pdf2[pI + 3] - feq);
542 408000 : boundary(_pdf2, pI, x, y, z, 3, _flag[index + nb], pI + 19 * nb);
543 408000 : feq -= 6.0 * _W[3] * _density[index] * cu;
544 408000 : _pdf2[pI + 15] -= _omega * (_pdf2[pI + 15] - feq);
545 408000 : boundary(_pdf2, pI, x, y, z, 15, _flag[index - nb], pI - 19 * nb);
546 : // pdf 4,14
547 408000 : cu = vel[1] - vel[2];
548 408000 : nb = _xO - _yO;
549 408000 : feq = _W[4] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
550 408000 : _pdf2[pI + 4] -= _omega * (_pdf2[pI + 4] - feq);
551 408000 : boundary(_pdf2, pI, x, y, z, 4, _flag[index + nb], pI + 19 * nb);
552 408000 : feq -= 6.0 * _W[4] * _density[index] * cu;
553 408000 : _pdf2[pI + 14] -= _omega * (_pdf2[pI + 14] - feq);
554 408000 : boundary(_pdf2, pI, x, y, z, 14, _flag[index - nb], pI - 19 * nb);
555 : // pdf 5,13
556 408000 : cu = -vel[0] - vel[1];
557 408000 : nb = -1 - _xO;
558 408000 : feq = _W[5] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
559 408000 : _pdf2[pI + 5] -= _omega * (_pdf2[pI + 5] - feq);
560 408000 : boundary(_pdf2, pI, x, y, z, 5, _flag[index + nb], pI + 19 * nb);
561 408000 : feq -= 6.0 * _W[5] * _density[index] * cu;
562 408000 : _pdf2[pI + 13] -= _omega * (_pdf2[pI + 13] - feq);
563 408000 : boundary(_pdf2, pI, x, y, z, 13, _flag[index - nb], pI - 19 * nb);
564 : // pdf 6,12
565 408000 : cu = -vel[1];
566 408000 : nb = -_xO;
567 408000 : feq = _W[6] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
568 408000 : _pdf2[pI + 6] -= _omega * (_pdf2[pI + 6] - feq);
569 408000 : boundary(_pdf2, pI, x, y, z, 6, _flag[index + nb], pI + 19 * nb);
570 408000 : feq -= 6.0 * _W[6] * _density[index] * cu;
571 408000 : _pdf2[pI + 12] -= _omega * (_pdf2[pI + 12] - feq);
572 408000 : boundary(_pdf2, pI, x, y, z, 12, _flag[index - nb], pI - 19 * nb);
573 : // pdf 7,11
574 408000 : cu = vel[0] - vel[1];
575 408000 : nb = 1 - _xO;
576 408000 : feq = _W[7] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
577 408000 : _pdf2[pI + 7] -= _omega * (_pdf2[pI + 7] - feq);
578 408000 : boundary(_pdf2, pI, x, y, z, 7, _flag[index + nb], pI + 19 * nb);
579 408000 : feq -= 6.0 * _W[7] * _density[index] * cu;
580 408000 : _pdf2[pI + 11] -= _omega * (_pdf2[pI + 11] - feq);
581 408000 : boundary(_pdf2, pI, x, y, z, 11, _flag[index - nb], pI - 19 * nb);
582 : // pdf 8,10
583 408000 : cu = -vel[0];
584 408000 : nb = -1;
585 408000 : feq = _W[8] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
586 408000 : _pdf2[pI + 8] -= _omega * (_pdf2[pI + 8] - feq);
587 408000 : boundary(_pdf2, pI, x, y, z, 8, _flag[index + nb], pI + 19 * nb);
588 408000 : feq -= 6.0 * _W[8] * _density[index] * cu;
589 408000 : _pdf2[pI + 10] -= _omega * (_pdf2[pI + 10] - feq);
590 408000 : boundary(_pdf2, pI, x, y, z, 10, _flag[index - nb], pI - 19 * nb);
591 : // pdf 9
592 408000 : _pdf2[pI + 9] -= _omega * (_pdf2[pI + 9] - _W[9] * _density[index] * u2);
593 408000 : }
594 :
595 : /** @brief takes care of the correct boundary treatment for the LB method
596 : * @param pdf particle distribution function
597 : * @param index start index for current cell in pdf-array
598 : * @param x the position in x direction of the cell
599 : * @param y the position in y direction of the cell
600 : * @param z the position in z direction of the cell
601 : * @param q distribution function number
602 : * @param flag boundary flag of neighbouring cell
603 : * @param nbIndex index of neighbouring cell */
604 7344000 : void boundary(double* const pdf, int index, int x, int y, int z, int q, const Flag& flag, int nbIndex) {
605 7344000 : if (flag != FLUID) {
606 599760 : if (flag == NO_SLIP) {
607 : // half-way bounce back
608 102000 : pdf[nbIndex + 18 - q] = pdf[index + q];
609 497760 : } else if (flag == MOVING_WALL) {
610 : // half-way bounce back + moving wall acceleration (only x-direction for
611 : // wall supported at the moment)
612 102000 : pdf[nbIndex + 18 - q] =
613 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]);
614 395760 : } else if (flag == PERIODIC) {
615 : // periodic treatment
616 0 : int target[3] = {x, y, z};
617 0 : if (target[0] + _C[q][0] == 0) {
618 0 : target[0] = _domainSizeX + 1;
619 0 : } else if (target[0] + _C[q][0] == _domainSizeX + 1) {
620 0 : target[0] = 0;
621 : }
622 0 : if (target[1] + _C[q][1] == 0) {
623 0 : target[1] = _domainSizeY + 1;
624 0 : } else if (target[1] + _C[q][1] == _domainSizeY + 1) {
625 0 : target[1] = 0;
626 : }
627 0 : if (target[2] + _C[q][2] == 0) {
628 0 : target[2] = _domainSizeZ + 1;
629 0 : } else if (target[2] + _C[q][2] == _domainSizeZ + 1) {
630 0 : target[2] = 0;
631 : }
632 0 : const int periodicNb = target[0] + (_domainSizeX + 2) * (target[1] + (_domainSizeY + 2) * target[2]);
633 0 : pdf[19 * periodicNb + q] = pdf[index + q];
634 : }
635 : }
636 7344000 : }
637 :
638 : /** @brief refers to the LB method; computes density and velocity on pdf
639 : * @param vel velocity
640 : * @param density density
641 : * @param pdf partial distribution function */
642 965240 : static void computeDensityAndVelocity(double* const vel, double& density, const double* const pdf) {
643 965240 : vel[0] = -(pdf[1] + pdf[5] + pdf[8] + pdf[11] + pdf[15]);
644 965240 : density = pdf[3] + pdf[7] + pdf[10] + pdf[13] + pdf[17];
645 965240 : vel[1] = (pdf[4] + pdf[11] + pdf[12] + pdf[13] + pdf[18]) - (pdf[0] + pdf[5] + pdf[6] + pdf[7] + pdf[14]);
646 965240 : vel[0] = density + vel[0];
647 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];
648 965240 : vel[2] = (pdf[14] + pdf[15] + pdf[16] + pdf[17] + pdf[18]) - (pdf[0] + pdf[1] + pdf[2] + pdf[3] + pdf[4]);
649 965240 : vel[0] = vel[0] / density;
650 965240 : vel[1] = vel[1] / density;
651 965240 : vel[2] = vel[2] / density;
652 965240 : }
653 :
654 : /** takes care of communication across one face in one direction.
655 : * @param pdf partial distribution function
656 : * @param sendBuffer send buffer
657 : * @param recvBuffer receive buffer
658 : * @param nbFlagTo direction into which message is sent
659 : * @param nbFlagFrom direction from which message is received
660 : * @param startSend 3d coordinates that define the start of the data to be
661 : * sent to neighbouring process
662 : * @param endSend 3d coordinates that define the end of the data to to be
663 : * sent to neighbouring process
664 : * @param startRecv 3d coordinates that define the start of the data to be
665 : * received from neighbouring process
666 : * @param endRecv 3d coordinates that define the end of the data to be
667 : * received from neighbouring process */
668 306 : void communicatePart(double* pdf, double* sendBuffer, double* recvBuffer, NbFlag nbFlagTo, NbFlag nbFlagFrom, tarch::la::Vector<3, int> startSend,
669 : tarch::la::Vector<3, int> endSend, tarch::la::Vector<3, int> startRecv, tarch::la::Vector<3, int> endRecv) {
670 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
671 : // directions that point to LEFT/RIGHT,... -> same ordering as enums!
672 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}};
673 306 : MPI_Request requests[2];
674 306 : MPI_Status status[2];
675 306 : tarch::la::Vector<2, int> plane;
676 306 : tarch::la::Vector<2, int> domainSize;
677 : // find out plane coordinates
678 306 : if (nbFlagTo == LEFT || nbFlagTo == RIGHT) {
679 102 : plane[0] = 1;
680 102 : plane[1] = 2;
681 102 : domainSize[0] = _domainSizeY;
682 102 : domainSize[1] = _domainSizeZ;
683 204 : } else if (nbFlagTo == FRONT || nbFlagTo == BACK) {
684 102 : plane[0] = 0;
685 102 : plane[1] = 2;
686 102 : domainSize[0] = _domainSizeX;
687 102 : domainSize[1] = _domainSizeZ;
688 102 : } else if (nbFlagTo == TOP || nbFlagTo == BOTTOM) {
689 102 : plane[0] = 0;
690 102 : plane[1] = 1;
691 102 : domainSize[0] = _domainSizeX;
692 102 : domainSize[1] = _domainSizeY;
693 : } else {
694 0 : std::cout << "ERROR LBCouetteSolver::communicatePart: d >2 or d < 0!" << std::endl;
695 0 : exit(EXIT_FAILURE);
696 : }
697 : // extract data and write to send buffer
698 306 : tarch::la::Vector<3, int> coords(0);
699 4488 : for (coords[2] = startSend[2]; coords[2] < endSend[2]; coords[2]++) {
700 49266 : for (coords[1] = startSend[1]; coords[1] < endSend[1]; coords[1]++) {
701 180132 : for (coords[0] = startSend[0]; coords[0] < endSend[0]; coords[0]++) {
702 810288 : for (int q = 0; q < 5; q++) {
703 675240 : sendBuffer[q + 5 * getParBuf(coords[plane[0]], coords[plane[1]], domainSize[0], domainSize[1])] =
704 675240 : pdf[directions[nbFlagTo][q] + 19 * get(coords[0], coords[1], coords[2])];
705 : }
706 : }
707 : }
708 : }
709 : // send and receive data
710 306 : MPI_Irecv(recvBuffer, (domainSize[0] + 2) * (domainSize[1] + 2) * 5, MPI_DOUBLE, _parallelNeighbours[nbFlagFrom], 1000,
711 306 : coupling::indexing::IndexingService<3>::getInstance().getComm(), &requests[0]);
712 306 : MPI_Isend(sendBuffer, (domainSize[0] + 2) * (domainSize[1] + 2) * 5, MPI_DOUBLE, _parallelNeighbours[nbFlagTo], 1000,
713 306 : coupling::indexing::IndexingService<3>::getInstance().getComm(), &requests[1]);
714 306 : MPI_Waitall(2, requests, status);
715 : // write data back to pdf field
716 306 : if (_parallelNeighbours[nbFlagFrom] != MPI_PROC_NULL) {
717 4284 : for (coords[2] = startRecv[2]; coords[2] < endRecv[2]; coords[2]++) {
718 46920 : for (coords[1] = startRecv[1]; coords[1] < endRecv[1]; coords[1]++) {
719 128520 : for (coords[0] = startRecv[0]; coords[0] < endRecv[0]; coords[0]++) {
720 514080 : for (int q = 0; q < 5; q++) {
721 428400 : if (_flag[get(coords[0], coords[1], coords[2])] == PARALLEL_BOUNDARY) {
722 428400 : pdf[directions[nbFlagTo][q] + 19 * get(coords[0], coords[1], coords[2])] =
723 428400 : recvBuffer[q + 5 * getParBuf(coords[plane[0]], coords[plane[1]], domainSize[0], domainSize[1])];
724 : }
725 : }
726 : }
727 : }
728 : }
729 : }
730 : #endif
731 306 : }
732 :
733 : /** @brief comunicates the boundary field data between the different processes
734 : */
735 51 : void communicate() {
736 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
737 : // send from right to left
738 102 : communicatePart(_pdf1, _sendBufferX, _recvBufferX, LEFT, RIGHT, tarch::la::Vector<3, int>(1, 1, 1),
739 51 : tarch::la::Vector<3, int>(2, _domainSizeY + 1, _domainSizeZ + 1), tarch::la::Vector<3, int>(_domainSizeX + 1, 1, 1),
740 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 1, _domainSizeZ + 1));
741 : // send from left to right
742 102 : communicatePart(_pdf1, _sendBufferX, _recvBufferX, RIGHT, LEFT, tarch::la::Vector<3, int>(_domainSizeX, 1, 1),
743 51 : tarch::la::Vector<3, int>(_domainSizeX + 1, _domainSizeY + 1, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, 1, 1),
744 51 : tarch::la::Vector<3, int>(1, _domainSizeY + 1, _domainSizeZ + 1));
745 : // send from back to front
746 102 : communicatePart(_pdf1, _sendBufferY, _recvBufferY, FRONT, BACK, tarch::la::Vector<3, int>(0, 1, 1),
747 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, 2, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, _domainSizeY + 1, 1),
748 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, _domainSizeZ + 1));
749 : // send from front to back
750 102 : communicatePart(_pdf1, _sendBufferY, _recvBufferY, BACK, FRONT, tarch::la::Vector<3, int>(0, _domainSizeY, 1),
751 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 1, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, 0, 1),
752 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, 1, _domainSizeZ + 1));
753 : // send from top to bottom
754 102 : communicatePart(_pdf1, _sendBufferZ, _recvBufferZ, BOTTOM, TOP, tarch::la::Vector<3, int>(0, 0, 1),
755 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, 2), tarch::la::Vector<3, int>(0, 0, _domainSizeZ + 1),
756 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, _domainSizeZ + 2));
757 : // send from bottom to top
758 102 : communicatePart(_pdf1, _sendBufferZ, _recvBufferZ, TOP, BOTTOM, tarch::la::Vector<3, int>(0, 0, _domainSizeZ),
759 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, 0, 0),
760 51 : tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, 1));
761 : #endif
762 51 : }
763 :
764 : /** @brief relaxation frequency */
765 : const double _omega;
766 : /** @brief velocity of moving wall of Couette flow */
767 : tarch::la::Vector<3, double> _wallVelocity;
768 : int _pdfsize{0};
769 : /** @brief partical distribution function field */
770 : double* _pdf1{NULL};
771 : /** @brief partial distribution function field (stores the old time step)*/
772 : double* _pdf2{NULL};
773 : /** @brief lattice velocities*/
774 : 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},
775 : {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}};
776 : /** @brief lattice weights */
777 : 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,
778 : 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};
779 : /** @brief enables avg_vel CSV output */
780 : const bool _plotAverageVelocity;
781 : };
782 :
783 : #endif // _MOLECULARDYNAMICS_COUPLING_SOLVERS_LBCOUETTESOLVER_H_
|
