/*  Lattice Boltzmann sample, written in C++, using the OpenLB
 *  library
 *
 *  Copyright (C) 2007, 2012 Jonas Latt, Mathias J. Krause
 *  Vojtech Cvrcek, Peter Weisbrod
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
 */

#include "olb2D.h"
#include "olb2D.hh"

#include <vector>

using namespace olb;
using namespace olb::descriptors;

typedef double T;

#define DESCRIPTOR D2Q9Descriptor

const T lx  = 4.;             // length of the channel
const T ly  = 1.;             // height of the channel
const int N = 30;             // resolution of the model
const T Re  = 100.;           // Reynolds number
const T maxPhysT = 60.;       // max. simulation time in s, SI unit
const T physInterval = 0.25;  // interval for the convergence check in s
const T residuum = 1e-5;      // residuum for the convergence check


void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
                     SuperGeometry2D<T>& superGeometry)
{

  OstreamManager clout(std::cout,"prepareGeometry");
  clout << "Prepare Geometry ..." << std::endl;

  superGeometry.rename(0,2);
  // Set material number for bounce back boundaries
  superGeometry.rename(2,1,0,1);

  T physSpacing = converter.getPhysDeltaX();
  Vector<T,2> extend{  physSpacing / 2, ly };
  Vector<T,2> origin{ -physSpacing / 4, 0  };

  // Set material number for inflow
  IndicatorCuboid2D<T> inflow(extend, origin);
  superGeometry.rename(1,3,1,inflow);

  // Set material number for outflow
  origin[0] = lx - physSpacing / 4;
  IndicatorCuboid2D<T> outflow(extend, origin);
  superGeometry.rename(1,4,1,outflow);

  IndicatorCuboid2D<T> obstacle(Vector<T,2> {0.2,0.1}, Vector<T,2> {1.25,0.45});
  superGeometry.rename(1,5,obstacle);

  // Removes all not needed boundary voxels outside the surface
  superGeometry.clean();
  // Removes all not needed boundary voxels inside the surface
  superGeometry.innerClean();
  superGeometry.checkForErrors();

  superGeometry.print();

  clout << "Prepare Geometry ... OK" << std::endl;
}

void prepareCoarseLattice(UnitConverter<T,DESCRIPTOR> const& converter,
                          SuperLattice2D<T, DESCRIPTOR>& sLattice,
                          Dynamics<T, DESCRIPTOR>& bulkDynamics,
                          sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition,
                          SuperGeometry2D<T>& superGeometry)
{
  OstreamManager clout(std::cout,"prepareLattice");
  clout << "Prepare coarse lattice ..." << std::endl;

  const T omega = converter.getLatticeRelaxationFrequency();

  sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
  sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk
  sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>());
  sLattice.defineDynamics(superGeometry, 3, &bulkDynamics); // inflow
  sLattice.defineDynamics(superGeometry, 4, &bulkDynamics); // outflow

  sBoundaryCondition.addVelocityBoundary(superGeometry, 3, omega);
  sBoundaryCondition.addPressureBoundary(superGeometry, 4, omega);

  const T Lx = converter.getLatticeLength(lx);
  const T Ly = converter.getLatticeLength(ly);
  const T p0 = 8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly);
  AnalyticalLinear2D<T,T> rho(-p0/lx*DESCRIPTOR<T>::invCs2, 0, p0*DESCRIPTOR<T>::invCs2+1);

  const T maxVelocity = converter.getCharLatticeVelocity();
  const T radius = ly/2;
  std::vector<T> axisPoint{0, ly/2};
  std::vector<T> axisDirection{1, 0};
  Poiseuille2D<T> u(axisPoint, axisDirection, maxVelocity, radius);

  sLattice.defineRhoU(superGeometry, 1, rho, u);
  sLattice.iniEquilibrium(superGeometry, 1, rho, u);
  sLattice.defineRhoU(superGeometry, 2, rho, u);
  sLattice.iniEquilibrium(superGeometry, 2, rho, u);

  sLattice.defineRhoU(superGeometry, 3, rho, u);
  sLattice.iniEquilibrium(superGeometry, 3, rho, u);
  sLattice.defineRhoU(superGeometry, 4, rho, u);
  sLattice.iniEquilibrium(superGeometry, 4, rho, u);

  sLattice.initialize();

  clout << "Prepare coarse lattice ... OK" << std::endl;
}

void prepareFineLattice(UnitConverter<T,DESCRIPTOR> const& converter,
                        SuperLattice2D<T, DESCRIPTOR>& sLattice,
                        Dynamics<T, DESCRIPTOR>& bulkDynamics,
                        sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition,
                        SuperGeometry2D<T>& superGeometry)
{
  OstreamManager clout(std::cout,"prepareLattice");
  clout << "Prepare fine lattice ..." << std::endl;

  const T omega = converter.getLatticeRelaxationFrequency();

  sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
  sLattice.defineDynamics(superGeometry, 1, &bulkDynamics); // bulk
  sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>());

  const T Lx = converter.getLatticeLength(lx);
  const T Ly = converter.getLatticeLength(ly);
  const T p0 = 8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly);
  AnalyticalLinear2D<T,T> rho(-p0/lx*DESCRIPTOR<T>::invCs2, 0, p0*DESCRIPTOR<T>::invCs2+1);

  const T maxVelocity = converter.getCharLatticeVelocity();
  const T radius = ly/2;
  std::vector<T> axisPoint{0, ly/2};
  std::vector<T> axisDirection{1, 0};
  Poiseuille2D<T> u(axisPoint, axisDirection, maxVelocity, radius);

  sLattice.defineRhoU(superGeometry, 1, rho, u);
  sLattice.iniEquilibrium(superGeometry, 1, rho, u);
  sLattice.defineRhoU(superGeometry, 2, rho, u);
  sLattice.iniEquilibrium(superGeometry, 2, rho, u);

  sLattice.initialize();

  clout << "Prepare fine lattice ... OK" << std::endl;
}

void getResults(const std::string& prefix,
                SuperLattice2D<T,DESCRIPTOR>& sLattice,
                UnitConverter<T,DESCRIPTOR> const& converter, int iT,
                SuperGeometry2D<T>& superGeometry, Timer<T>& timer, bool hasConverged)
{
  OstreamManager clout(std::cout,"getResults");

  SuperVTMwriter2D<T> vtmWriter(prefix + "poiseuille2d");
  SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity(sLattice, converter);
  SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure(sLattice, converter);
  SuperLatticeGeometry2D<T, DESCRIPTOR> geometry( sLattice, superGeometry );
  vtmWriter.addFunctor(geometry);
  vtmWriter.addFunctor(velocity);
  vtmWriter.addFunctor(pressure);

  const int statIter = converter.getLatticeTime(maxPhysT/10.);

  if (iT==0) {
    vtmWriter.createMasterFile();
  }

  if (iT%100==0) {
    vtmWriter.write(iT);
  }

  if (iT%statIter==0 || hasConverged) {
    timer.update(iT);
    timer.printStep();

    sLattice.getStatistics().print(iT,converter.getPhysTime(iT));
  }
}

template <typename T, template<typename> class DESCRIPTOR> class FineCoupler;
template <typename T, template<typename> class DESCRIPTOR> class CoarseCoupler;

template <typename T, template<typename> class DESCRIPTOR>
class Grid {
private:
  std::unique_ptr<UnitConverter<T,DESCRIPTOR>>  _converter;
  std::unique_ptr<CuboidGeometry2D<T>>          _cuboids;
  std::unique_ptr<LoadBalancer<T>>              _balancer;
  std::unique_ptr<SuperGeometry2D<T>>           _geometry;
  std::unique_ptr<SuperLattice2D<T,DESCRIPTOR>> _lattice;

  std::vector<std::unique_ptr<Grid<T,DESCRIPTOR>>> _fineGrids;

  std::vector<std::unique_ptr<FineCoupler<T,DESCRIPTOR>>> _fineCouplers;
  std::vector<std::unique_ptr<CoarseCoupler<T,DESCRIPTOR>>> _coarseCouplers;

public:
  static std::unique_ptr<Grid<T,DESCRIPTOR>> make(IndicatorF2D<T>& domainF, int resolution, T tau, int re)
  {
    return std::unique_ptr<Grid<T,DESCRIPTOR>>(
             new Grid<T,DESCRIPTOR>(domainF, resolution, tau, re)
           );
  }

  Grid(IndicatorF2D<T>& domainF, int resolution, T tau, int re):
    _converter(new UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR>(
                 resolution, // resolution: number of voxels per charPhysL
                 tau,        // latticeRelaxationTime: relaxation time, has to be greater than 0.5!
                 T{1},       // charPhysLength: reference length of simulation geometry
                 T{1},       // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
                 T{1./re},   // physViscosity: physical kinematic viscosity in __m^2 / s__
                 T{1},       // physDensity: physical density in __kg / m^3__
                 T{1})),
    _cuboids(new CuboidGeometry2D<T>(
               domainF,
               _converter->getConversionFactorLength(),
               1)),
    _balancer(new HeuristicLoadBalancer<T>(
                *_cuboids)),
    _geometry(new SuperGeometry2D<T>(
                *_cuboids,
                *_balancer,
                2)),
    _lattice(new SuperLattice2D<T,DESCRIPTOR>(
               *_geometry))
  {
    _converter->print();
  }

  UnitConverter<T,DESCRIPTOR>& getConverter()
  {
    return *_converter;
  }

  CuboidGeometry2D<T>& getCuboidGeometry()
  {
    return *_cuboids;
  }

  LoadBalancer<T>& getLoadBalancer()
  {
    return *_balancer;
  }

  SuperGeometry2D<T>& getSuperGeometry()
  {
    return *_geometry;
  }

  SuperLattice2D<T,DESCRIPTOR>& getSuperLattice()
  {
    return *_lattice;
  }

  T getScalingFactor()
  {
    return 4. - _converter->getLatticeRelaxationFrequency();
  }

  T getInvScalingFactor()
  {
    return 1./getScalingFactor();
  }

  void collideAndStream()
  {
    for ( auto& fineCoupler : _fineCouplers ) {
      fineCoupler->store();
    }

    this->getSuperLattice().collideAndStream();

    for ( auto& fineGrid : _fineGrids ) {
      fineGrid->getSuperLattice().collideAndStream();
    }

    for ( auto& fineCoupler : _fineCouplers ) {
      fineCoupler->interpolate();
      fineCoupler->couple();
    }

    for ( auto& fineGrid : _fineGrids ) {
      fineGrid->getSuperLattice().collideAndStream();
    }

    for ( auto& fineCoupler : _fineCouplers ) {
      fineCoupler->store();
      fineCoupler->couple();
    }

    for ( auto& coarseCoupler : _coarseCouplers ) {
      coarseCoupler->couple();
    }
  }

  Grid<T,DESCRIPTOR>* refine(Vector<T,2> origin, Vector<T,2> extend)
  {
    IndicatorCuboid2D<T> fineCuboid(extend, origin);
    _fineGrids.emplace_back(
      new Grid<T,DESCRIPTOR>(
        fineCuboid,
        2*getConverter().getResolution(),
        2.0*getConverter().getLatticeRelaxationTime() - 0.5,
        getConverter().getReynoldsNumber()
      ));
    Grid<T,DESCRIPTOR>* const fineGrid = _fineGrids.back().get();

    const T coarseDeltaX = getConverter().getPhysDeltaX();

    _fineCouplers.emplace_back(
      new FineCoupler<T,DESCRIPTOR>(
        *this, *fineGrid, origin, Vector<T,2> {0,extend[1]}));
    _fineCouplers.emplace_back(
      new FineCoupler<T,DESCRIPTOR>(
        *this, *fineGrid, origin + Vector<T,2> {extend[0],0}, Vector<T,2> {0,extend[1]}));
    _fineCouplers.emplace_back(
      new FineCoupler<T,DESCRIPTOR>(
        *this, *fineGrid, origin + Vector<T,2>(0,extend[1]), Vector<T,2> {extend[0],0}));
    _fineCouplers.emplace_back(
      new FineCoupler<T,DESCRIPTOR>(
        *this, *fineGrid, origin, Vector<T,2> {extend[0],0}));

    _coarseCouplers.emplace_back(
      new CoarseCoupler<T,DESCRIPTOR>(
        *this, *fineGrid, origin + coarseDeltaX, Vector<T,2> {0,extend[1]-2*coarseDeltaX}));
    _coarseCouplers.emplace_back(
      new CoarseCoupler<T,DESCRIPTOR>(
        *this, *fineGrid, origin + Vector<T,2>(extend[0]-coarseDeltaX,coarseDeltaX), Vector<T,2> {0,extend[1]-2*coarseDeltaX}));
    _coarseCouplers.emplace_back(
      new CoarseCoupler<T,DESCRIPTOR>(
        *this, *fineGrid, origin + Vector<T,2> {coarseDeltaX,extend[1]-coarseDeltaX}, Vector<T,2> {extend[0]-2*coarseDeltaX,0}));
    _coarseCouplers.emplace_back(
      new CoarseCoupler<T,DESCRIPTOR>(
        *this, *fineGrid, origin + coarseDeltaX, Vector<T,2> {extend[0]-2*coarseDeltaX,0}));

    return fineGrid;
  }
};

template <typename T, template<typename> class DESCRIPTOR>
void computeFeq(const Cell<T,DESCRIPTOR>& cell, T fEq[DESCRIPTOR<T>::q])
{
  T rho{};
  T u[2] {};
  cell.computeRhoU(rho, u);
  const T uSqr = u[0]*u[0] + u[1]*u[1];
  for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
    fEq[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr);
  }
}

template <typename T, template<typename> class DESCRIPTOR>
void computeFneq(const Cell<T,DESCRIPTOR>& cell, T fNeq[DESCRIPTOR<T>::q])
{
  T rho{};
  T u[2] {};
  cell.computeRhoU(rho, u);
  lbHelpers<T,DESCRIPTOR>::computeFneq(cell, fNeq, rho, u);
}

T order4interpolation(T fm3, T fm1, T f1, T f3)
{
  return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3);
}

T order3interpolation(T fm1, T f1, T f3)
{
  return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3;
}

T order2interpolation(T f0, T f1)
{
  return 0.5 * (f0 + f1);
}

template <typename T, template<typename> class DESCRIPTOR>
class Coupler {
protected:
  Grid<T,DESCRIPTOR>& _coarse;
  Grid<T,DESCRIPTOR>& _fine;

  const int  _coarseSize;
  const int  _fineSize;
  const bool _vertical;

  Vector<T,2>   _physOrigin;
  Vector<int,3> _coarseOrigin;
  Vector<int,3> _fineOrigin;

  Vector<int,3> getFineLatticeR(int y) const
  {
    if (_vertical) {
      return _fineOrigin + Vector<int,3> {0, 0, y};
    }
    else {
      return _fineOrigin + Vector<int,3> {0, y, 0};
    }
  }

  Vector<int,3> getCoarseLatticeR(int y) const
  {
    if (_vertical) {
      return _coarseOrigin + Vector<int,3> {0, 0, y};
    }
    else {
      return _coarseOrigin + Vector<int,3> {0, y, 0};
    }
  }

public:
  Coupler(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine, Vector<T,2> origin, Vector<T,2> extend):
    _coarse(coarse),
    _fine(fine),
    _coarseSize(floor(extend.norm() / coarse.getConverter().getPhysDeltaX() + 0.5)+1),
    _fineSize(2*_coarseSize-1),
    _vertical(util::nearZero(extend[0]))
  {
    OLB_ASSERT(util::nearZero(extend[0]) || util::nearZero(extend[1]), "Coupling is only implemented alongside unit vectors");

    const auto& coarseGeometry = _coarse.getCuboidGeometry();
    const auto& fineGeometry   = _fine.getCuboidGeometry();

    Vector<int,3> tmpLatticeOrigin{};
    coarseGeometry.getLatticeR(origin, tmpLatticeOrigin);
    _physOrigin = coarseGeometry.getPhysR(tmpLatticeOrigin.toStdVector());

    coarseGeometry.getLatticeR(_physOrigin, _coarseOrigin);
    fineGeometry.getLatticeR(_physOrigin, _fineOrigin);
  }
};

template <typename T, template<typename> class DESCRIPTOR>
class FineCoupler : public Coupler<T,DESCRIPTOR> {
private:
  std::vector<T>           _c2f_rho;
  std::vector<Vector<T,2>> _c2f_u;
  std::vector<Vector<T,DESCRIPTOR<T>::q>> _c2f_fneq;

public:
  FineCoupler(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine, Vector<T,2> origin, Vector<T,2> extend):
    Coupler<T,DESCRIPTOR>(coarse, fine, origin, extend),
    _c2f_rho(this->_coarseSize),
    _c2f_u(this->_coarseSize, Vector<T,2>(T{})),
    _c2f_fneq(this->_coarseSize, Vector<T,DESCRIPTOR<T>::q>(T{}))
  {
    OstreamManager clout(std::cout,"C2F");
    clout << "coarse origin: " << this->_coarseOrigin[0] << " " << this->_coarseOrigin[1] << " " << this->_coarseOrigin[2] << std::endl;
    clout << "fine origin:   " << this->_fineOrigin[0] << " " << this->_fineOrigin[1] << " " << this->_fineOrigin[2] << std::endl;
    clout << "fine size:   " << this->_fineSize << std::endl;
  }

  void store()
  {
    auto& coarseLattice = this->_coarse.getSuperLattice();

    for (int y=0; y < this->_coarseSize; ++y) {
      const auto pos = this->getCoarseLatticeR(y);

      T rho{};
      T u[2] {};
      T fNeq[DESCRIPTOR<T>::q] {};
      coarseLattice.get(pos).computeRhoU(rho, u);
      computeFneq(coarseLattice.get(pos), fNeq);

      _c2f_rho[y]  = rho;
      _c2f_u[y]    = Vector<T,2>(u);
      _c2f_fneq[y] = Vector<T,DESCRIPTOR<T>::q>(fNeq);
    }
  }

  void interpolate()
  {
    auto& coarseLattice = this->_coarse.getSuperLattice();

    for (int y=0; y < this->_coarseSize; ++y) {
      const auto coarsePos = this->getCoarseLatticeR(y);

      T rho{};
      T u[2] {};
      coarseLattice.get(coarsePos).computeRhoU(rho, u);
      _c2f_rho[y]  = order2interpolation(rho, _c2f_rho[y]);
      _c2f_u[y][0] = order2interpolation(u[0], _c2f_u[y][0]);
      _c2f_u[y][1] = order2interpolation(u[1], _c2f_u[y][1]);

      T fNeq[DESCRIPTOR<T>::q] {};
      computeFneq(coarseLattice.get(coarsePos), fNeq);
      for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
        _c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], _c2f_fneq[y][iPop]);
      }
    }
  }

  void couple()
  {
    auto& coarseLattice = this->_coarse.getSuperLattice();
    auto& fineLattice   = this->_fine.getSuperLattice();

    for (int y=0; y < this->_coarseSize; ++y) {
      const auto coarsePos = this->getCoarseLatticeR(y);
      const auto finePos   = this->getFineLatticeR(2*y);

      T fEq[DESCRIPTOR<T>::q] {};
      computeFeq(coarseLattice.get(coarsePos), fEq);

      for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
        fineLattice.get(finePos)[iPop] = fEq[iPop] + this->_coarse.getScalingFactor() * _c2f_fneq[y][iPop];
      }
    }

    for (int y=1; y < this->_coarseSize-2; ++y) {
      const T rho = order4interpolation(
                      _c2f_rho[y-1],
                      _c2f_rho[y],
                      _c2f_rho[y+1],
                      _c2f_rho[y+2]
                    );
      T u[2] {};
      u[0] = order4interpolation(
               _c2f_u[y-1][0],
               _c2f_u[y][0],
               _c2f_u[y+1][0],
               _c2f_u[y+2][0]
             );
      u[1] = order4interpolation(
               _c2f_u[y-1][1],
               _c2f_u[y][1],
               _c2f_u[y+1][1],
               _c2f_u[y+2][1]
             );
      const T uSqr = u[0]*u[0] + u[1]*u[1];

      const auto finePos = this->getFineLatticeR(1+2*y);

      for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
        const T fneq = order4interpolation(
                         _c2f_fneq[y-1][iPop],
                         _c2f_fneq[y][iPop],
                         _c2f_fneq[y+1][iPop],
                         _c2f_fneq[y+2][iPop]
                       );

        fineLattice.get(finePos)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
      }
    }

    {
      const T rho = order3interpolation(
                      _c2f_rho[0],
                      _c2f_rho[1],
                      _c2f_rho[2]
                    );
      T u[2] {};
      u[0] = order3interpolation(
               _c2f_u[0][0],
               _c2f_u[1][0],
               _c2f_u[2][0]
             );
      u[1] = order3interpolation(
               _c2f_u[0][1],
               _c2f_u[1][1],
               _c2f_u[2][1]
             );
      const T uSqr = u[0]*u[0] + u[1]*u[1];

      const auto finePos = this->getFineLatticeR(1);

      for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
        const T fneq = order3interpolation(
                         _c2f_fneq[0][iPop],
                         _c2f_fneq[1][iPop],
                         _c2f_fneq[2][iPop]
                       );
        fineLattice.get(finePos)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
      }
    }

    {
      const T rho = order3interpolation(
                      _c2f_rho[this->_coarseSize-1],
                      _c2f_rho[this->_coarseSize-2],
                      _c2f_rho[this->_coarseSize-3]
                    );
      T u[2] {};
      u[0] = order3interpolation(
               _c2f_u[this->_coarseSize-1][0],
               _c2f_u[this->_coarseSize-2][0],
               _c2f_u[this->_coarseSize-3][0]
             );
      u[1] = order3interpolation(
               _c2f_u[this->_coarseSize-1][1],
               _c2f_u[this->_coarseSize-2][1],
               _c2f_u[this->_coarseSize-3][1]
             );
      const T uSqr = u[0]*u[0] + u[1]*u[1];

      const auto finePos = this->getFineLatticeR(this->_fineSize-2);

      for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
        const T fneq = order3interpolation(
                         _c2f_fneq[this->_coarseSize-1][iPop],
                         _c2f_fneq[this->_coarseSize-2][iPop],
                         _c2f_fneq[this->_coarseSize-3][iPop]
                       );
        fineLattice.get(finePos)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
      }
    }
  }
};

template <typename T, template<typename> class DESCRIPTOR>
void computeRestrictedFneq(const SuperLattice2D<T,DESCRIPTOR>& lattice, Vector<int,3> pos, T restrictedFneq[DESCRIPTOR<T>::q])
{
  for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
    T fNeq[DESCRIPTOR<T>::q] {};
    computeFneq(lattice.get(pos[0], pos[1] + DESCRIPTOR<T>::c[iPop][0], pos[2] + DESCRIPTOR<T>::c[iPop][1]), fNeq);

    for (int jPop=0; jPop < DESCRIPTOR<T>::q; ++jPop) {
      restrictedFneq[jPop] += fNeq[jPop];
    }
  }

  for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
    restrictedFneq[iPop] /= DESCRIPTOR<T>::q;
  }
}

template <typename T, template<typename> class DESCRIPTOR>
class CoarseCoupler : public Coupler<T,DESCRIPTOR> {
public:
  CoarseCoupler(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine, Vector<T,2> origin, Vector<T,2> extend):
    Coupler<T,DESCRIPTOR>(coarse, fine, origin, extend)
  {
    OstreamManager clout(std::cout,"F2C");
    clout << "coarse origin: " << this->_coarseOrigin[0] << " " << this->_coarseOrigin[1] << " " << this->_coarseOrigin[2] << std::endl;
    clout << "fine origin:   " << this->_fineOrigin[0] << " " << this->_fineOrigin[1] << " " << this->_fineOrigin[2] << std::endl;
    clout << "coarse size:   " << this->_coarseSize << std::endl;
  }

  void couple()
  {
    auto& fineLattice = this->_fine.getSuperLattice();
    auto& coarseLattice = this->_coarse.getSuperLattice();

    for (int y=0; y < this->_coarseSize; ++y) {
      const auto finePos   = this->getFineLatticeR(2*y);
      const auto coarsePos = this->getCoarseLatticeR(y);

      T fEq[DESCRIPTOR<T>::q] {};
      computeFeq(fineLattice.get(finePos), fEq);

      T fNeq[DESCRIPTOR<T>::q] {};
      computeRestrictedFneq(fineLattice, finePos, fNeq);

      for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
        coarseLattice.get(coarsePos)[iPop] = fEq[iPop] + this->_coarse.getInvScalingFactor() * fNeq[iPop];
      }
    }
  }
};

int main(int argc, char* argv[])
{
  olbInit(&argc, &argv);
  singleton::directories().setOutputDir("./tmp/");
  OstreamManager clout(std::cout,"main");

  const T coarseDeltaX = 1./N;

  Vector<T,2> coarseOrigin { 0.0, 0.0 };
  Vector<T,2> coarseExtend { lx, ly  };
  IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);

  Vector<T,2> fineOrigin   { 0.25 * lx, 0.2 };
  Vector<T,2> fineExtend   { 0.6  * lx, 0.6 };

  auto coarseGrid = Grid<T,DESCRIPTOR>::make(coarseCuboid, N, 0.8, Re);
  auto fineGrid   = coarseGrid->refine(fineOrigin, fineExtend);

  prepareGeometry(coarseGrid->getConverter(), coarseGrid->getSuperGeometry());
  prepareGeometry(fineGrid->getConverter(), fineGrid->getSuperGeometry());

  fineGrid->getSuperGeometry().rename(2,1);
  fineGrid->getSuperGeometry().rename(5,2);
  IndicatorCuboid2D<T> overlapCuboid(fineExtend - 4*coarseDeltaX, fineOrigin + 2*coarseDeltaX);
  coarseGrid->getSuperGeometry().rename(1,0,overlapCuboid);
  coarseGrid->getSuperGeometry().rename(2,0,overlapCuboid);
  coarseGrid->getSuperGeometry().rename(5,0,overlapCuboid);

  Dynamics<T, DESCRIPTOR>* coarseBulkDynamics;
  coarseBulkDynamics = new BGKdynamics<T, DESCRIPTOR>(
    coarseGrid->getConverter().getLatticeRelaxationFrequency(),
    instances::getBulkMomenta<T, DESCRIPTOR>());

  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> coarseSBoundaryCondition(coarseGrid->getSuperLattice());
  createLocalBoundaryCondition2D<T, DESCRIPTOR>(coarseSBoundaryCondition);

  prepareCoarseLattice(
    coarseGrid->getConverter(),
    coarseGrid->getSuperLattice(),
    *coarseBulkDynamics,
    coarseSBoundaryCondition,
    coarseGrid->getSuperGeometry());

  Dynamics<T, DESCRIPTOR>* fineBulkDynamics;
  fineBulkDynamics = new BGKdynamics<T, DESCRIPTOR>(
    fineGrid->getConverter().getLatticeRelaxationFrequency(),
    instances::getBulkMomenta<T, DESCRIPTOR>());

  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid->getSuperLattice());
  createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition);

  prepareFineLattice(
    fineGrid->getConverter(),
    fineGrid->getSuperLattice(),
    *fineBulkDynamics,
    fineSBoundaryCondition,
    fineGrid->getSuperGeometry());

  clout << "starting simulation..." << endl;
  Timer<T> timer(
    coarseGrid->getConverter().getLatticeTime(maxPhysT),
    coarseGrid->getSuperGeometry().getStatistics().getNvoxel());
  util::ValueTracer<T> converge(
    fineGrid->getConverter().getLatticeTime(physInterval),
    residuum);
  timer.start();

  for (int iT = 0; iT < coarseGrid->getConverter().getLatticeTime(maxPhysT); ++iT) {
    if (converge.hasConverged()) {
      clout << "Simulation converged." << endl;
      break;
    }

    coarseGrid->collideAndStream();

    getResults(
      "coarse_",
      coarseGrid->getSuperLattice(),
      coarseGrid->getConverter(),
      iT,
      coarseGrid->getSuperGeometry(),
      timer,
      converge.hasConverged());
    getResults(
      "fine_",
      fineGrid->getSuperLattice(),
      fineGrid->getConverter(),
      iT,
      fineGrid->getSuperGeometry(),
      timer,
      converge.hasConverged());

    converge.takeValue(fineGrid->getSuperLattice().getStatistics().getAverageEnergy(), true);
  }

  timer.stop();
  timer.printSummary();
}