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