Add hacky MPI support for grid refinement

Works but is nowhere near anything one could consider good.
Obvious issues:

* More than one cuboid per grid makes it harder to determine the next lattice cell to be coupled
	* i.e. currently lattice positions are determined ad hoc by resolving their physical position
* Coupling is not actually parallelized
	* All coupling lines are traversed by all processes, way to much communication
* Load balancing and cuboid decomposition doesn't care about refinement
	* ideally refined cuboids should be computationally near their coarse _parent_ cuboids

The first two isses should be fixable with a reasonable amount of work.
This sadly doesn't apply in any form to the last issue.

2 files changed, 65 insertions, 19 deletions

diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh
index 8370707..ab082f4 100644
--- a/src/refinement/coupler2D.hh
+++ b/src/refinement/coupler2D.hh
@@ -35,10 +35,16 @@ template <typename T, template<typename> class DESCRIPTOR>
 Vector<int,3> Coupler2D<T,DESCRIPTOR>::getFineLatticeR(int y) const
 {
   if (_vertical) {
-    return _fineOrigin + Vector<int,3> {0, 0, y};
+    const auto physR = _physOrigin + Vector<T,2>{0, y*_fine.getConverter().getPhysDeltaX()};
+    Vector<int,3> latticeR{};
+    _fine.getCuboidGeometry().getLatticeR(physR, latticeR);
+    return latticeR;
   }
   else {
-    return _fineOrigin + Vector<int,3> {0, y, 0};
+    const auto physR = _physOrigin + Vector<T,2>{y*_fine.getConverter().getPhysDeltaX(), 0};
+    Vector<int,3> latticeR{};
+    _fine.getCuboidGeometry().getLatticeR(physR, latticeR);
+    return latticeR;
   }
 }
 
@@ -46,10 +52,16 @@ template <typename T, template<typename> class DESCRIPTOR>
 Vector<int,3> Coupler2D<T,DESCRIPTOR>::getCoarseLatticeR(int y) const
 {
   if (_vertical) {
-    return _coarseOrigin + Vector<int,3> {0, 0, y};
+    const auto physR = _physOrigin + Vector<T,2>{0, y*_coarse.getConverter().getPhysDeltaX()};
+    Vector<int,3> latticeR{};
+    _coarse.getCuboidGeometry().getLatticeR(physR, latticeR);
+    return latticeR;
   }
   else {
-    return _coarseOrigin + Vector<int,3> {0, y, 0};
+    const auto physR = _physOrigin + Vector<T,2>{y*_coarse.getConverter().getPhysDeltaX(), 0};
+    Vector<int,3> latticeR{};
+    _coarse.getCuboidGeometry().getLatticeR(physR, latticeR);
+    return latticeR;
   }
 }
 
@@ -97,12 +109,13 @@ void FineCoupler2D<T,DESCRIPTOR>::store()
 
   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);
-    coarseLattice.get(pos).computeFneq(fNeq);
+    Cell<T,DESCRIPTOR> coarseCell;
+    coarseLattice.get(pos, coarseCell);
+    lbHelpers<T,DESCRIPTOR>::computeRhoU(coarseCell, rho, u);
+    lbHelpers<T,DESCRIPTOR>::computeFneq(coarseCell, fNeq, rho, u);
 
     _c2f_rho[y]  = rho;
     _c2f_u[y]    = Vector<T,2>(u);
@@ -134,17 +147,18 @@ void FineCoupler2D<T,DESCRIPTOR>::interpolate()
   auto& coarseLattice = this->_coarse.getSuperLattice();
 
   for (int y=0; y < this->_coarseSize; ++y) {
-    const auto coarseCell = coarseLattice.get(this->getCoarseLatticeR(y));
+    Cell<T,DESCRIPTOR> coarseCell;
+    coarseLattice.get(this->getCoarseLatticeR(y), coarseCell);
 
     T rho{};
     T u[2] {};
-    coarseCell.computeRhoU(rho, u);
+    lbHelpers<T,DESCRIPTOR>::computeRhoU(coarseCell, 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] {};
-    coarseCell.computeFneq(fNeq);
+    lbHelpers<T,DESCRIPTOR>::computeFneq(coarseCell, fNeq, rho, u);
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
       _c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], _c2f_fneq[y][iPop]);
     }
@@ -162,11 +176,16 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
     const auto finePos   = this->getFineLatticeR(2*y);
 
     T fEq[DESCRIPTOR<T>::q] {};
-    coarseLattice.get(coarsePos).computeFeq(fEq);
+    Cell<T,DESCRIPTOR> coarseCell;
+    coarseLattice.get(coarsePos, coarseCell);
+    lbHelpers<T,DESCRIPTOR>::computeFeq(coarseCell, fEq);
 
+    Cell<T,DESCRIPTOR> cell;
+    fineLattice.get(finePos, cell);
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      fineLattice.get(finePos)[iPop] = fEq[iPop] + this->_coarse.getScalingFactor() * _c2f_fneq[y][iPop];
+       cell[iPop] = fEq[iPop] + this->_coarse.getScalingFactor() * _c2f_fneq[y][iPop];
     }
+    fineLattice.set(finePos, cell);
   }
 
   for (int y=1; y < this->_coarseSize-2; ++y) {
@@ -192,6 +211,8 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
     const T uSqr = u[0]*u[0] + u[1]*u[1];
 
     const auto finePos = this->getFineLatticeR(1+2*y);
+    Cell<T,DESCRIPTOR> fineCell;
+    fineLattice.get(finePos, fineCell);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
       const T fneq = order4interpolation(
@@ -201,8 +222,10 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
                        _c2f_fneq[y+2][iPop]
                      );
 
-      fineLattice.get(finePos)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
+      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
     }
+
+    fineLattice.set(finePos, fineCell);
   }
 
   {
@@ -225,6 +248,8 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
     const T uSqr = u[0]*u[0] + u[1]*u[1];
 
     const auto finePos = this->getFineLatticeR(1);
+    Cell<T,DESCRIPTOR> fineCell;
+    fineLattice.get(finePos, fineCell);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
       const T fneq = order3interpolation(
@@ -232,8 +257,10 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
                        _c2f_fneq[1][iPop],
                        _c2f_fneq[2][iPop]
                      );
-      fineLattice.get(finePos)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
+      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
     }
+
+    fineLattice.set(finePos, fineCell);
   }
 
   {
@@ -256,6 +283,8 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
     const T uSqr = u[0]*u[0] + u[1]*u[1];
 
     const auto finePos = this->getFineLatticeR(this->_fineSize-2);
+    Cell<T,DESCRIPTOR> fineCell;
+    fineLattice.get(finePos, fineCell);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
       const T fneq = order3interpolation(
@@ -263,8 +292,10 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
                        _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;
+      fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
     }
+
+    fineLattice.set(finePos, fineCell);
   }
 }
 
@@ -276,7 +307,10 @@ void computeRestrictedFneq(const SuperLattice2D<T,DESCRIPTOR>& lattice,
 {
   for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
     T fNeq[DESCRIPTOR<T>::q] {};
-    lattice.get(pos[0], pos[1] + DESCRIPTOR<T>::c[iPop][0], pos[2] + DESCRIPTOR<T>::c[iPop][1]).computeFneq(fNeq);
+    Cell<T,DESCRIPTOR> cell;
+    const Vector<int,3> neighborPos {pos[0], pos[1] + DESCRIPTOR<T>::c[iPop][0], pos[2] + DESCRIPTOR<T>::c[iPop][1]};
+    lattice.get(neighborPos, cell);
+    lbHelpers<T,DESCRIPTOR>::computeFneq(cell, fNeq);
 
     for (int jPop=0; jPop < DESCRIPTOR<T>::q; ++jPop) {
       restrictedFneq[jPop] += fNeq[jPop];
@@ -311,14 +345,21 @@ void CoarseCoupler2D<T,DESCRIPTOR>::couple()
     const auto coarsePos = this->getCoarseLatticeR(y);
 
     T fEq[DESCRIPTOR<T>::q] {};
-    fineLattice.get(finePos).computeFeq(fEq);
+    Cell<T,DESCRIPTOR> fineCell;
+    fineLattice.get(finePos, fineCell);
+    lbHelpers<T,DESCRIPTOR>::computeFeq(fineCell, fEq);
 
     T fNeq[DESCRIPTOR<T>::q] {};
     computeRestrictedFneq(fineLattice, finePos, fNeq);
 
+    Cell<T,DESCRIPTOR> coarseCell;
+    coarseLattice.get(coarsePos, coarseCell);
+
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      coarseLattice.get(coarsePos)[iPop] = fEq[iPop] + this->_coarse.getInvScalingFactor() * fNeq[iPop];
+      coarseCell[iPop] = fEq[iPop] + this->_coarse.getInvScalingFactor() * fNeq[iPop];
     }
+
+    coarseLattice.set(coarsePos, coarseCell);
   }
 }
 
diff --git a/src/refinement/grid2D.hh b/src/refinement/grid2D.hh
index 71b943e..65dd84d 100644
--- a/src/refinement/grid2D.hh
+++ b/src/refinement/grid2D.hh
@@ -55,7 +55,12 @@ Grid2D<T,DESCRIPTOR>::Grid2D(IndicatorF2D<T>& domainF, int resolution, T tau, in
   _cuboids(new CuboidGeometry2D<T>(
              domainF,
              _converter->getConversionFactorLength(),
-             1)),
+#ifdef PARALLEL_MODE_MPI
+             singleton::mpi().getSize()
+#else
+             1
+#endif
+             )),
   _balancer(new HeuristicLoadBalancer<T>(
               *_cuboids)),
   _geometry(new SuperGeometry2D<T>(