LBCouetteSolver.h

// Copyright (C) 2015 Technische Universitaet Muenchen
// This file is part of the Mamico project. For conditions of distribution
// and use, please see the copyright notice in Mamico's main folder
#ifndef _MOLECULARDYNAMICS_COUPLING_SOLVERS_LBCOUETTESOLVER_H_
#define _MOLECULARDYNAMICS_COUPLING_SOLVERS_LBCOUETTESOLVER_H_
 
namespace coupling {
namespace solvers {
class LBCouetteSolver;
class LBCouetteSolverState;
} // namespace solvers
} // namespace coupling
 
#if defined(_OPENMP)
#include <omp.h>
#endif
#include "coupling/interface/PintableMacroSolver.h"
#include "coupling/solvers/NumericalSolver.h"
#include <cmath>
 
class coupling::solvers::LBCouetteSolverState : public coupling::interface::PintableMacroSolverState {
public:
  LBCouetteSolverState(int size) : _pdf(size, 0) {}
 
  LBCouetteSolverState(int size, double* pdf) : LBCouetteSolverState(size) { std::copy(pdf, pdf + size, _pdf.data()); }
 
  std::unique_ptr<State> clone() const override { return std::make_unique<LBCouetteSolverState>(*this); }
 
  ~LBCouetteSolverState() {}
 
  int getSizeBytes() const override { return sizeof(double) * _pdf.size(); }
 
  std::unique_ptr<State> operator+(const State& rhs) override;
  std::unique_ptr<State> operator-(const State& rhs) override;
 
  bool operator==(const State& rhs) const override {
    const LBCouetteSolverState* other = dynamic_cast<const LBCouetteSolverState*>(&rhs);
    if (other == nullptr)
      return false;
    return _pdf == other->_pdf;
  }
 
  double* getData() override { return _pdf.data(); }
  const double* getData() const override { return _pdf.data(); }
 
  void print(std::ostream& os) const override { os << "<LBCouetteSolverState instance with size " << getSizeBytes() << ">"; }
 
private:
  std::vector<double> _pdf;
};

class coupling::solvers::LBCouetteSolver : public coupling::solvers::NumericalSolver, public coupling::interface::PintableMacroSolver {
public:
  LBCouetteSolver(const double channelheight, tarch::la::Vector<3, double> wallVelocity, const double kinVisc, const double dx, const double dt,
                  const int plotEveryTimestep, const bool plotAverageVelocity, const std::string filestem, const tarch::la::Vector<3, unsigned int> processes,
                  const unsigned int numThreads = 1, const Scenario* scen = nullptr)
      : coupling::solvers::NumericalSolver(channelheight, dx, dt, kinVisc, plotEveryTimestep, filestem, processes, scen), _mode(Mode::coupling), _dt_pint(dt),
        _omega(1.0 / (3.0 * (kinVisc * dt / (dx * dx)) + 0.5)), _wallVelocity((dt / dx) * wallVelocity), _plotAverageVelocity(plotAverageVelocity) {
    // return if required
    if (skipRank()) {
      return;
    }
    _pdfsize = 19 * (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2);
    _pdf1 = new double[_pdfsize];
    _pdf2 = new double[_pdfsize];
#if defined(_OPENMP)
    omp_set_num_threads(numThreads);
#endif
#if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
    std::cout << "Domain size=" << _domainSizeX << "," << _domainSizeY << "," << _domainSizeZ << std::endl;
    std::cout << "tau=" << 1.0 / _omega << std::endl;
    std::cout << "wallVelocity=" << _wallVelocity << std::endl;
    for (int z = 0; z < _domainSizeZ + 2; z++) {
      for (int y = 0; y < _domainSizeY + 2; y++) {
        for (int x = 0; x < _domainSizeX + 2; x++) {
          std::cout << x << "," << y << "," << z << "FLAG=" << _flag[get(x, y, z)] << std::endl;
        }
      }
    }
#endif
    // check pointers
    if ((_pdf1 == NULL) || (_pdf2 == NULL) || (_vel == NULL) || (_density == NULL) || (_flag == NULL)) {
      std::cout << "ERROR LBCouetteSolver: NULL ptr!" << std::endl;
      exit(EXIT_FAILURE);
    }
#if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
    if ((_sendBufferX == NULL) || (_recvBufferX == NULL) || (_sendBufferY == NULL) || (_recvBufferY == NULL) || (_sendBufferZ == NULL) ||
        (_recvBufferZ == NULL)) {
      std::cout << "ERROR LBCouetteSolver: NULL ptr in send/recv!" << std::endl;
      exit(EXIT_FAILURE);
    }
#endif
// init everything with lattice weights
#pragma omp parallel for
    for (int i = 0; i < (_domainSizeX + 2) * (_domainSizeY + 2) * (_domainSizeZ + 2); i++) {
      for (int q = 0; q < 19; q++) {
        _pdf1[get(i) * 19 + q] = _W[q];
        _pdf2[get(i) * 19 + q] = _W[q];
      }
    }
    computeDensityAndVelocityEverywhere();
  }

  virtual ~LBCouetteSolver() {
    if (_pdf1 != NULL) {
      delete[] _pdf1;
      _pdf1 = NULL;
    }
    if (_pdf2 != NULL) {
      delete[] _pdf2;
      _pdf2 = NULL;
    }
    if (_vel != NULL) {
      delete[] _vel;
      _vel = NULL;
    }
    if (_density != NULL) {
      delete[] _density;
      _density = NULL;
    }
    if (_flag != NULL) {
      delete[] _flag;
      _flag = NULL;
    }
#if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
    if (_sendBufferX != NULL) {
      delete[] _sendBufferX;
      _sendBufferX = NULL;
    }
    if (_sendBufferY != NULL) {
      delete[] _sendBufferY;
      _sendBufferY = NULL;
    }
    if (_sendBufferZ != NULL) {
      delete[] _sendBufferZ;
      _sendBufferZ = NULL;
    }
    if (_recvBufferX != NULL) {
      delete[] _recvBufferX;
      _recvBufferX = NULL;
    }
    if (_recvBufferY != NULL) {
      delete[] _recvBufferY;
      _recvBufferY = NULL;
    }
    if (_recvBufferZ != NULL) {
      delete[] _recvBufferZ;
      _recvBufferZ = NULL;
    }
#endif
  }

  void advance(double dt) override {
    if (skipRank()) {
      return;
    }
    const int timesteps = floor(dt / _dt + 0.5);
    if (fabs(timesteps * _dt - dt) / _dt > 1.0e-8) {
      std::cout << "ERROR LBCouetteSolver::advance(): time steps and dt do not match!" << std::endl;
      exit(EXIT_FAILURE);
    }
    for (int i = 0; i < timesteps; i++) {
      if (_plotEveryTimestep >= 1 && _counter % _plotEveryTimestep == 0)
        computeDensityAndVelocityEverywhere();
      plot();
      plot_avg_vel();
      collidestream();
      communicate(); // exchange between neighbouring MPI subdomains
      _counter++;
    }
  }

  void setMDBoundaryValues(coupling::datastructures::FlexibleCellContainer<3>& md2macroBuffer) override {
    if (skipRank()) {
      return;
    }
#if (COUPLING_MD_ERROR == COUPLING_MD_YES)
    if (_mode == Mode::supervising) {
      std::cout << "ERROR LBCouetteSolver setMDBoundaryValues() called in supervising mode" << std::endl;
      exit(EXIT_FAILURE);
    }
#endif
 
    // loop over all received cells
    for (auto pair : md2macroBuffer) {
      I01 idx;
      const coupling::datastructures::CouplingCell<3>* couplingCell;
      std::tie(couplingCell, idx) = pair;
      // determine cell index of this cell in LB domain
      tarch::la::Vector<3, unsigned int> globalCellCoords{idx.get()};
      globalCellCoords[0] = (globalCellCoords[0] + _offset[0]) - _coords[0] * _avgDomainSizeX;
      globalCellCoords[1] = (globalCellCoords[1] + _offset[1]) - _coords[1] * _avgDomainSizeY;
      globalCellCoords[2] = (globalCellCoords[2] + _offset[2]) - _coords[2] * _avgDomainSizeZ;
#if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
      std::cout << "Process coords: " << _coords << ":  GlobalCellCoords for index " << idx << ": " << globalCellCoords << std::endl;
#endif
      const int index = get(globalCellCoords[0], globalCellCoords[1], globalCellCoords[2]);
#if (COUPLING_MD_ERROR == COUPLING_MD_YES)
      if (_flag[index] != MD_BOUNDARY) {
        std::cout << "ERROR LBCouetteSolver::setMDBoundaryValues(): Cell " << index << " is no MD boundary cell!" << std::endl;
        exit(EXIT_FAILURE);
      }
#endif
      // set velocity value and pdfs in MD boundary cell (before streaming); the
      // boundary velocities are interpolated between the neighbouring and this
      // cell. This interpolation is valid for FLUID-MD_BOUNDARY neighbouring
      // relations only. determine local velocity received from MaMiCo and
      // convert it to LB units; store the velocity in _vel
      tarch::la::Vector<3, double> localVel((1.0 / couplingCell->getMacroscopicMass()) * (_dt / _dx) * couplingCell->getMacroscopicMomentum());
      for (unsigned int d = 0; d < 3; d++) {
        _vel[3 * index + d] = localVel[d];
      }
      // loop over all pdfs and set them according to interpolated moving-wall
      // conditions
      for (unsigned int q = 0; q < 19; q++) {
        // index of neighbour cell; only if cell is located inside local domain
        if (((int)globalCellCoords[0] + _C[q][0] > 0) && ((int)globalCellCoords[0] + _C[q][0] < _domainSizeX + 1) &&
            ((int)globalCellCoords[1] + _C[q][1] > 0) && ((int)globalCellCoords[1] + _C[q][1] < _domainSizeY + 1) &&
            ((int)globalCellCoords[2] + _C[q][2] > 0) && ((int)globalCellCoords[2] + _C[q][2] < _domainSizeZ + 1)) {
          const int nbIndex = get((_C[q][0] + globalCellCoords[0]), (_C[q][1] + globalCellCoords[1]), (_C[q][2] + globalCellCoords[2]));
          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]),
                                                         0.5 * (_vel[3 * index + 2] + _vel[3 * nbIndex + 2]));
          _pdf1[19 * index + q] =
              _pdf1[19 * nbIndex + 18 - q] -
              6.0 * _W[q] * _density[nbIndex] * (_C[18 - q][0] * interpolVel[0] + _C[18 - q][1] * interpolVel[1] + _C[18 - q][2] * interpolVel[2]);
        }
      }
    }
  }

  tarch::la::Vector<3, double> getVelocity(tarch::la::Vector<3, double> pos) const override {
    tarch::la::Vector<3, unsigned int> coords;
    const tarch::la::Vector<3, double> domainOffset(_coords[0] * _dx * _avgDomainSizeX, _coords[1] * _dx * _avgDomainSizeY, _coords[2] * _dx * _avgDomainSizeZ);
    // check pos-data for process locality (todo: put this in debug mode in
    // future releases)
    if ((pos[0] < domainOffset[0]) || (pos[0] > domainOffset[0] + _domainSizeX * _dx) || (pos[1] < domainOffset[1]) ||
        (pos[1] > domainOffset[1] + _domainSizeY * _dx) || (pos[2] < domainOffset[2]) || (pos[2] > domainOffset[2] + _domainSizeZ * _dx)) {
      std::cout << "ERROR LBCouetteSolver::getVelocity(): Position " << pos << " out of range!" << std::endl;
      std::cout << "domainOffset = " << domainOffset << std::endl;
      std::cout << "_domainSizeX = " << _domainSizeX << std::endl;
      std::cout << "_domainSizeY = " << _domainSizeY << std::endl;
      std::cout << "_domainSizeZ = " << _domainSizeZ << std::endl;
      std::cout << "_dx = " << _dx << std::endl;
      exit(EXIT_FAILURE);
    }
    // compute index for respective cell (_dx+... for ghost cells); use coords
    // to store local cell coordinates
    for (unsigned int d = 0; d < 3; d++) {
      coords[d] = (unsigned int)((_dx + pos[d] - domainOffset[d]) / _dx);
    }
    const int index = get(coords[0], coords[1], coords[2]);
    tarch::la::Vector<3, double> vel(0.0);
    // extract and scale velocity to "real"=MD units
    for (int d = 0; d < 3; d++) {
      vel[d] = _dx / _dt * _vel[3 * index + d];
    }
#if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
    std::cout << "Position " << pos << " corresponds to cell: " << coords << "; vel=" << vel << std::endl;
#endif
    return vel;
  }

  double getDensity(tarch::la::Vector<3, double> pos) const override {
    tarch::la::Vector<3, unsigned int> coords;
    const tarch::la::Vector<3, double> domainOffset(_coords[0] * _dx * _avgDomainSizeX, _coords[1] * _dx * _avgDomainSizeY, _coords[2] * _dx * _avgDomainSizeZ);
    // check pos-data for process locality (todo: put this in debug mode in
    // future releases)
    if ((pos[0] < domainOffset[0]) || (pos[0] > domainOffset[0] + _domainSizeX * _dx) || (pos[1] < domainOffset[1]) ||
        (pos[1] > domainOffset[1] + _domainSizeY * _dx) || (pos[2] < domainOffset[2]) || (pos[2] > domainOffset[2] + _domainSizeZ * _dx)) {
      std::cout << "ERROR LBCouetteSolver::getDensity(): Position " << pos << " out of range!" << std::endl;
      exit(EXIT_FAILURE);
    }
    // compute index for respective cell (_dx+... for ghost cells); use coords
    // to store local cell coordinates
    for (unsigned int d = 0; d < 3; d++) {
      coords[d] = (unsigned int)((_dx + pos[d] - domainOffset[d]) / _dx);
    }
    const int index = get(coords[0], coords[1], coords[2]);
    return _density[index];
  }

  virtual void setWallVelocity(const tarch::la::Vector<3, double> wallVelocity) override { _wallVelocity = (_dt / _dx) * wallVelocity; }

 
  std::unique_ptr<State> getState() override {
    computeDensityAndVelocityEverywhere();
    return std::make_unique<LBCouetteSolverState>(_pdfsize, _pdf1);
  }
 
  void setState(const std::unique_ptr<State>& input, int cycle) override {
    const LBCouetteSolverState* state = dynamic_cast<const LBCouetteSolverState*>(input.get());
 
#if (COUPLING_MD_ERROR == COUPLING_MD_YES)
    if (state == nullptr) {
      std::cout << "ERROR LBCouetteSolver setState() wrong state type" << std::endl;
      exit(EXIT_FAILURE);
    }
#endif
 
    std::copy(state->getData(), state->getData() + _pdfsize, _pdf1);
    computeDensityAndVelocityEverywhere();
 
    _counter = cycle;
  }
 
  std::unique_ptr<State> operator()(const std::unique_ptr<State>& input, int cycle) override {
    setState(input, cycle);
 
#if (COUPLING_MD_ERROR == COUPLING_MD_YES)
    if (_mode != Mode::supervising) {
      std::cout << "ERROR LBCouetteSolver operator() called but not in supervising mode" << std::endl;
      exit(EXIT_FAILURE);
    }
#endif
 
    advance(_dt_pint);
    return getState();
  }
 
  Mode getMode() const override { return _mode; }

  std::unique_ptr<PintableMacroSolver> getSupervisor(int num_cycles, double visc_multiplier) const override {
#if (COUPLING_MD_ERROR == COUPLING_MD_YES)
    if (_mode == Mode::supervising) {
      std::cout << "ERROR LBCouetteSolver getSupervisor(): already in supervising mode" << std::endl;
      exit(EXIT_FAILURE);
    }
#endif
 
    int numThreads = 1;
#if defined(_OPENMP)
    numThreads = omp_get_num_threads();
#endif
 
    auto res = std::make_unique<LBCouetteSolver>(_channelheight, _wallVelocity * _dx / _dt, _kinVisc * visc_multiplier, _dx, _dt, _plotEveryTimestep, _plotAverageVelocity,
                                                 _filestem + std::string("_supervising"), _processes, numThreads, _scen);
 
    res->_mode = Mode::supervising;
    res->_dt_pint = _dt * num_cycles;
 
    return res;
  }
 
  void print(std::ostream& os) const override {
    if (_mode == Mode::supervising)
      os << "<LBCouetteSolver instance in supervising mode >";
    if (_mode == Mode::coupling)
      os << "<LBCouetteSolver instance in coupling mode >";
  }
 
  double get_avg_vel(const std::unique_ptr<State>& state) const override {
    double vel[3];
    double density;
    double res[3]{0, 0, 0};
    for (int i = 0; i < _pdfsize; i += 19) {
      LBCouetteSolver::computeDensityAndVelocity(vel, density, state->getData() + i);
      res[0] += vel[0];
      res[1] += vel[1];
      res[2] += vel[2];
    }
    if (_pdfsize > 0) {
      res[0] /= (_pdfsize / 19);
      res[1] /= (_pdfsize / 19);
      res[2] /= (_pdfsize / 19);
    }
    return std::sqrt(res[0] * res[0] + res[1] * res[1] + res[2] * res[2]);
  }
 
private:
  Mode _mode;
  double _dt_pint;
 
  void computeDensityAndVelocityEverywhere() {
    if (skipRank())
      return;
    for (int z = 1; z < _domainSizeZ + 1; z++) {
      for (int y = 1; y < _domainSizeY + 1; y++) {
        for (int x = 1; x < _domainSizeX + 1; x++) {
          const int index = get(x, y, z);
          const int pI = 19 * index;
          double* vel = &_vel[3 * index];
          computeDensityAndVelocity(vel, _density[index], &_pdf1[pI]);
        }
      }
    }
  }
 
  void plot_avg_vel(){
    if (!_plotAverageVelocity) return;
 
    int rank = 0;
#if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
    MPI_Comm_rank(coupling::indexing::IndexingService<3>::getInstance().getComm(), &rank);
#endif
    std::stringstream ss;
    ss << _filestem << "_r" << rank;
    if (_scen != nullptr) {
      auto ts = _scen->getTimeIntegrationService();
      if (ts != nullptr) {
        if (ts->isPintEnabled())
          ss << "_i" << ts->getInteration();
      }
    }
    ss << ".csv";
    std::string filename = ss.str();
    std::ofstream file(filename.c_str(), _counter==0 ? std::ofstream::out : std::ofstream::app);
    if (!file.is_open()) {
      std::cout << "ERROR LBCouetteSolver::plot_avg_vel(): Could not open file " << filename << "!" << std::endl;
      exit(EXIT_FAILURE);
    }
 
    if(_counter == 0){
      file << "coupling_cycle ; avg_vel" << std::endl;
    }
 
    std::unique_ptr<State> s = std::make_unique<LBCouetteSolverState>(_pdfsize, _pdf1);
    double vel = get_avg_vel(s);
    file << _counter << " ; " << vel << std::endl;
    file.close();
  }


  void collidestream() {
#pragma omp parallel for
    for (int z = 1; z < _domainSizeZ + 1; z++) {
      for (int y = 1; y < _domainSizeY + 1; y++) {
        for (int x = 1; x < _domainSizeX + 1; x++) {
          const int index = get(x, y, z);
          if (_flag[index] == FLUID) {
            stream(index);
            collide(index, x, y, z);
          }
        }
      }
    }
    // swap fields
    double* swap = _pdf1;
    _pdf1 = _pdf2;
    _pdf2 = swap;
  }

  void stream(int index) {
    const int pI = 19 * index;
    for (int q = 0; q < 9; q++) {
      const int nb = 19 * (_C[q][0] + _C[q][1] * _xO + _C[q][2] * _yO);
      _pdf2[pI + q] = _pdf1[pI + q - nb];
      _pdf2[pI + 18 - q] = _pdf1[pI + 18 - q + nb];
    }
    _pdf2[pI + 9] = _pdf1[pI + 9];
  }

  void collide(int index, int x, int y, int z) {
    // index of start of cell-local pdfs in AoS
    const int pI = 19 * index;
    // compute and store density, velocity
    double* vel = &_vel[3 * index];
    computeDensityAndVelocity(vel, _density[index], &_pdf2[pI]);
    // collide (BGK); always handle pdfs no. q and inv(q)=18-q in one step
    const double u2 = 1.0 - 1.5 * (vel[0] * vel[0] + vel[1] * vel[1] + vel[2] * vel[2]);
    // pdf 0,18
    double cu = -vel[1] - vel[2];
    int nb = -_xO - _yO;
    double feq = _W[0] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
    _pdf2[pI] -= _omega * (_pdf2[pI] - feq);
    boundary(_pdf2, pI, x, y, z, 0, _flag[index + nb], pI + 19 * nb);
    feq -= 6.0 * _W[0] * _density[index] * cu;
    _pdf2[pI + 18] -= _omega * (_pdf2[pI + 18] - feq);
    boundary(_pdf2, pI, x, y, z, 18, _flag[index - nb], pI - 19 * nb);
    // pdf 1,17
    cu = -vel[0] - vel[2];
    nb = -1 - _yO;
    feq = _W[1] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
    _pdf2[pI + 1] -= _omega * (_pdf2[pI + 1] - feq);
    boundary(_pdf2, pI, x, y, z, 1, _flag[index + nb], pI + 19 * nb);
    feq -= 6.0 * _W[1] * _density[index] * cu;
    _pdf2[pI + 17] -= _omega * (_pdf2[pI + 17] - feq);
    boundary(_pdf2, pI, x, y, z, 17, _flag[index - nb], pI - 19 * nb);
    // pdf 2,16
    cu = -vel[2];
    nb = -_yO;
    feq = _W[2] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
    _pdf2[pI + 2] -= _omega * (_pdf2[pI + 2] - feq);
    boundary(_pdf2, pI, x, y, z, 2, _flag[index + nb], pI + 19 * nb);
    feq -= 6.0 * _W[2] * _density[index] * cu;
    _pdf2[pI + 16] -= _omega * (_pdf2[pI + 16] - feq);
    boundary(_pdf2, pI, x, y, z, 16, _flag[index - nb], pI - 19 * nb);
    // pdf 3,15
    cu = vel[0] - vel[2];
    nb = 1 - _yO;
    feq = _W[3] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
    _pdf2[pI + 3] -= _omega * (_pdf2[pI + 3] - feq);
    boundary(_pdf2, pI, x, y, z, 3, _flag[index + nb], pI + 19 * nb);
    feq -= 6.0 * _W[3] * _density[index] * cu;
    _pdf2[pI + 15] -= _omega * (_pdf2[pI + 15] - feq);
    boundary(_pdf2, pI, x, y, z, 15, _flag[index - nb], pI - 19 * nb);
    // pdf 4,14
    cu = vel[1] - vel[2];
    nb = _xO - _yO;
    feq = _W[4] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
    _pdf2[pI + 4] -= _omega * (_pdf2[pI + 4] - feq);
    boundary(_pdf2, pI, x, y, z, 4, _flag[index + nb], pI + 19 * nb);
    feq -= 6.0 * _W[4] * _density[index] * cu;
    _pdf2[pI + 14] -= _omega * (_pdf2[pI + 14] - feq);
    boundary(_pdf2, pI, x, y, z, 14, _flag[index - nb], pI - 19 * nb);
    // pdf 5,13
    cu = -vel[0] - vel[1];
    nb = -1 - _xO;
    feq = _W[5] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
    _pdf2[pI + 5] -= _omega * (_pdf2[pI + 5] - feq);
    boundary(_pdf2, pI, x, y, z, 5, _flag[index + nb], pI + 19 * nb);
    feq -= 6.0 * _W[5] * _density[index] * cu;
    _pdf2[pI + 13] -= _omega * (_pdf2[pI + 13] - feq);
    boundary(_pdf2, pI, x, y, z, 13, _flag[index - nb], pI - 19 * nb);
    // pdf 6,12
    cu = -vel[1];
    nb = -_xO;
    feq = _W[6] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
    _pdf2[pI + 6] -= _omega * (_pdf2[pI + 6] - feq);
    boundary(_pdf2, pI, x, y, z, 6, _flag[index + nb], pI + 19 * nb);
    feq -= 6.0 * _W[6] * _density[index] * cu;
    _pdf2[pI + 12] -= _omega * (_pdf2[pI + 12] - feq);
    boundary(_pdf2, pI, x, y, z, 12, _flag[index - nb], pI - 19 * nb);
    // pdf 7,11
    cu = vel[0] - vel[1];
    nb = 1 - _xO;
    feq = _W[7] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
    _pdf2[pI + 7] -= _omega * (_pdf2[pI + 7] - feq);
    boundary(_pdf2, pI, x, y, z, 7, _flag[index + nb], pI + 19 * nb);
    feq -= 6.0 * _W[7] * _density[index] * cu;
    _pdf2[pI + 11] -= _omega * (_pdf2[pI + 11] - feq);
    boundary(_pdf2, pI, x, y, z, 11, _flag[index - nb], pI - 19 * nb);
    // pdf 8,10
    cu = -vel[0];
    nb = -1;
    feq = _W[8] * _density[index] * (u2 + 3.0 * cu + 4.5 * cu * cu);
    _pdf2[pI + 8] -= _omega * (_pdf2[pI + 8] - feq);
    boundary(_pdf2, pI, x, y, z, 8, _flag[index + nb], pI + 19 * nb);
    feq -= 6.0 * _W[8] * _density[index] * cu;
    _pdf2[pI + 10] -= _omega * (_pdf2[pI + 10] - feq);
    boundary(_pdf2, pI, x, y, z, 10, _flag[index - nb], pI - 19 * nb);
    // pdf 9
    _pdf2[pI + 9] -= _omega * (_pdf2[pI + 9] - _W[9] * _density[index] * u2);
  }

  void boundary(double* const pdf, int index, int x, int y, int z, int q, const Flag& flag, int nbIndex) {
    if (flag != FLUID) {
      if (flag == NO_SLIP) {
        // half-way bounce back
        pdf[nbIndex + 18 - q] = pdf[index + q];
      } else if (flag == MOVING_WALL) {
        // half-way bounce back + moving wall acceleration (only x-direction for
        // wall supported at the moment)
        pdf[nbIndex + 18 - q] =
            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]);
      } else if (flag == PERIODIC) {
        // periodic treatment
        int target[3] = {x, y, z};
        if (target[0] + _C[q][0] == 0) {
          target[0] = _domainSizeX + 1;
        } else if (target[0] + _C[q][0] == _domainSizeX + 1) {
          target[0] = 0;
        }
        if (target[1] + _C[q][1] == 0) {
          target[1] = _domainSizeY + 1;
        } else if (target[1] + _C[q][1] == _domainSizeY + 1) {
          target[1] = 0;
        }
        if (target[2] + _C[q][2] == 0) {
          target[2] = _domainSizeZ + 1;
        } else if (target[2] + _C[q][2] == _domainSizeZ + 1) {
          target[2] = 0;
        }
        const int periodicNb = target[0] + (_domainSizeX + 2) * (target[1] + (_domainSizeY + 2) * target[2]);
        pdf[19 * periodicNb + q] = pdf[index + q];
      }
    }
  }

  static void computeDensityAndVelocity(double* const vel, double& density, const double* const pdf) {
    vel[0] = -(pdf[1] + pdf[5] + pdf[8] + pdf[11] + pdf[15]);
    density = pdf[3] + pdf[7] + pdf[10] + pdf[13] + pdf[17];
    vel[1] = (pdf[4] + pdf[11] + pdf[12] + pdf[13] + pdf[18]) - (pdf[0] + pdf[5] + pdf[6] + pdf[7] + pdf[14]);
    vel[0] = density + vel[0];
    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];
    vel[2] = (pdf[14] + pdf[15] + pdf[16] + pdf[17] + pdf[18]) - (pdf[0] + pdf[1] + pdf[2] + pdf[3] + pdf[4]);
    vel[0] = vel[0] / density;
    vel[1] = vel[1] / density;
    vel[2] = vel[2] / density;
  }

  void communicatePart(double* pdf, double* sendBuffer, double* recvBuffer, NbFlag nbFlagTo, NbFlag nbFlagFrom, tarch::la::Vector<3, int> startSend,
                       tarch::la::Vector<3, int> endSend, tarch::la::Vector<3, int> startRecv, tarch::la::Vector<3, int> endRecv) {
#if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
    // directions that point to LEFT/RIGHT,... -> same ordering as enums!
    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}};
    MPI_Request requests[2];
    MPI_Status status[2];
    tarch::la::Vector<2, int> plane;
    tarch::la::Vector<2, int> domainSize;
    // find out plane coordinates
    if (nbFlagTo == LEFT || nbFlagTo == RIGHT) {
      plane[0] = 1;
      plane[1] = 2;
      domainSize[0] = _domainSizeY;
      domainSize[1] = _domainSizeZ;
    } else if (nbFlagTo == FRONT || nbFlagTo == BACK) {
      plane[0] = 0;
      plane[1] = 2;
      domainSize[0] = _domainSizeX;
      domainSize[1] = _domainSizeZ;
    } else if (nbFlagTo == TOP || nbFlagTo == BOTTOM) {
      plane[0] = 0;
      plane[1] = 1;
      domainSize[0] = _domainSizeX;
      domainSize[1] = _domainSizeY;
    } else {
      std::cout << "ERROR LBCouetteSolver::communicatePart: d >2 or d < 0!" << std::endl;
      exit(EXIT_FAILURE);
    }
    // extract data and write to send buffer
    tarch::la::Vector<3, int> coords(0);
    for (coords[2] = startSend[2]; coords[2] < endSend[2]; coords[2]++) {
      for (coords[1] = startSend[1]; coords[1] < endSend[1]; coords[1]++) {
        for (coords[0] = startSend[0]; coords[0] < endSend[0]; coords[0]++) {
          for (int q = 0; q < 5; q++) {
            sendBuffer[q + 5 * getParBuf(coords[plane[0]], coords[plane[1]], domainSize[0], domainSize[1])] =
                pdf[directions[nbFlagTo][q] + 19 * get(coords[0], coords[1], coords[2])];
          }
        }
      }
    }
    // send and receive data
    MPI_Irecv(recvBuffer, (domainSize[0] + 2) * (domainSize[1] + 2) * 5, MPI_DOUBLE, _parallelNeighbours[nbFlagFrom], 1000,
              coupling::indexing::IndexingService<3>::getInstance().getComm(), &requests[0]);
    MPI_Isend(sendBuffer, (domainSize[0] + 2) * (domainSize[1] + 2) * 5, MPI_DOUBLE, _parallelNeighbours[nbFlagTo], 1000,
              coupling::indexing::IndexingService<3>::getInstance().getComm(), &requests[1]);
    MPI_Waitall(2, requests, status);
    // write data back to pdf field
    if (_parallelNeighbours[nbFlagFrom] != MPI_PROC_NULL) {
      for (coords[2] = startRecv[2]; coords[2] < endRecv[2]; coords[2]++) {
        for (coords[1] = startRecv[1]; coords[1] < endRecv[1]; coords[1]++) {
          for (coords[0] = startRecv[0]; coords[0] < endRecv[0]; coords[0]++) {
            for (int q = 0; q < 5; q++) {
              if (_flag[get(coords[0], coords[1], coords[2])] == PARALLEL_BOUNDARY) {
                pdf[directions[nbFlagTo][q] + 19 * get(coords[0], coords[1], coords[2])] =
                    recvBuffer[q + 5 * getParBuf(coords[plane[0]], coords[plane[1]], domainSize[0], domainSize[1])];
              }
            }
          }
        }
      }
    }
#endif
  }

  void communicate() {
#if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
    // send from right to left
    communicatePart(_pdf1, _sendBufferX, _recvBufferX, LEFT, RIGHT, tarch::la::Vector<3, int>(1, 1, 1),
                    tarch::la::Vector<3, int>(2, _domainSizeY + 1, _domainSizeZ + 1), tarch::la::Vector<3, int>(_domainSizeX + 1, 1, 1),
                    tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 1, _domainSizeZ + 1));
    // send from left to right
    communicatePart(_pdf1, _sendBufferX, _recvBufferX, RIGHT, LEFT, tarch::la::Vector<3, int>(_domainSizeX, 1, 1),
                    tarch::la::Vector<3, int>(_domainSizeX + 1, _domainSizeY + 1, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, 1, 1),
                    tarch::la::Vector<3, int>(1, _domainSizeY + 1, _domainSizeZ + 1));
    // send from back to front
    communicatePart(_pdf1, _sendBufferY, _recvBufferY, FRONT, BACK, tarch::la::Vector<3, int>(0, 1, 1),
                    tarch::la::Vector<3, int>(_domainSizeX + 2, 2, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, _domainSizeY + 1, 1),
                    tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, _domainSizeZ + 1));
    // send from front to back
    communicatePart(_pdf1, _sendBufferY, _recvBufferY, BACK, FRONT, tarch::la::Vector<3, int>(0, _domainSizeY, 1),
                    tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 1, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, 0, 1),
                    tarch::la::Vector<3, int>(_domainSizeX + 2, 1, _domainSizeZ + 1));
    // send from top to bottom
    communicatePart(_pdf1, _sendBufferZ, _recvBufferZ, BOTTOM, TOP, tarch::la::Vector<3, int>(0, 0, 1),
                    tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, 2), tarch::la::Vector<3, int>(0, 0, _domainSizeZ + 1),
                    tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, _domainSizeZ + 2));
    // send from bottom to top
    communicatePart(_pdf1, _sendBufferZ, _recvBufferZ, TOP, BOTTOM, tarch::la::Vector<3, int>(0, 0, _domainSizeZ),
                    tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, _domainSizeZ + 1), tarch::la::Vector<3, int>(0, 0, 0),
                    tarch::la::Vector<3, int>(_domainSizeX + 2, _domainSizeY + 2, 1));
#endif
  }

  const double _omega;
  tarch::la::Vector<3, double> _wallVelocity;
  int _pdfsize{0};
  double* _pdf1{NULL};
  double* _pdf2{NULL};
  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},
                      {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}};

  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,
                      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};

  const bool _plotAverageVelocity;
};
 
#endif // _MOLECULARDYNAMICS_COUPLING_SOLVERS_LBCOUETTESOLVER_H_