Refactor grid coupling into classes

1 files changed, 245 insertions, 172 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 22c35cb..632fba5 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -38,7 +38,7 @@ typedef double T;
 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 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
@@ -325,182 +325,237 @@ T order2interpolation(T f0, T f1)
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
-void storeC2F(
-  Grid<T,DESCRIPTOR>& coarse,
-  T c2f_rho[],
-  T c2f_u[][2],
-  T c2f_fneq[][DESCRIPTOR<T>::q]
-)
-{
-  auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
-
-  const int x = coarseLattice.getNx()-2;
-  for (int y=0; y < coarseLattice.getNy(); ++y) {
-    T rho{};
-    T u[2] {};
-    coarseLattice.get(x,y).computeRhoU(rho, u);
-    c2f_rho[y]  = rho;
-    c2f_u[y][0] = u[0];
-    c2f_u[y][1] = u[1];
-
-    T fNeq[DESCRIPTOR<T>::q] {};
-    computeFneq(coarseLattice.get(x,y), fNeq);
-    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      c2f_fneq[y][iPop] = fNeq[iPop];
+class Coupler {
+protected:
+  Grid<T,DESCRIPTOR>& _coarse;
+  Grid<T,DESCRIPTOR>& _fine;
+  const Vector<T,2> _origin;
+  const Vector<T,2> _extend;
+  const int _coarseSize;
+  const bool _vertical;
+  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};
     }
   }
-}
 
-template <typename T, template<typename> class DESCRIPTOR>
-void interpolateStoredC2F(
-  Grid<T,DESCRIPTOR>& coarse,
-  T c2f_rho[],
-  T c2f_u[][2],
-  T c2f_fneq[][DESCRIPTOR<T>::q]
-)
-{
-  auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
-
-  const int x = coarseLattice.getNx()-2;
-  for (int y=0; y < coarseLattice.getNy(); ++y) {
-    T rho{};
-    T u[2] {};
-    coarseLattice.get(x,y).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(x,y), fNeq);
-    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], c2f_fneq[y][iPop]);
+  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};
     }
   }
-}
 
-template <typename T, template<typename> class DESCRIPTOR>
-void coupleC2F(
-  Grid<T,DESCRIPTOR>& coarse,
-  Grid<T,DESCRIPTOR>& fine,
-  const T c2f_rho[],
-  const T c2f_u[][2],
-  const T c2f_fneq[][DESCRIPTOR<T>::q]
-)
-{
-  auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
-  auto& fineLattice   = fine.getSuperLattice().getBlockLattice(0);
+public:
+  Coupler(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine, Vector<T,2> origin, Vector<T,2> extend):
+    _coarse(coarse),
+    _fine(fine),
+    _origin(origin),
+    _extend(extend),
+    _coarseSize(extend.norm() / coarse.getConverter().getPhysDeltaX()),
+    _vertical(util::nearZero(extend[0]))
+  {
+    OLB_ASSERT(util::nearZero(extend[0]) || util::nearZero(extend[1]), "Coupling is only implemented alongside unit vectors");
 
-  const int x = coarseLattice.getNx()-2;
+    _coarse.getCuboidGeometry().getFloorLatticeR(_origin, _coarseOrigin);
+    _fine.getCuboidGeometry().getFloorLatticeR(_origin, _fineOrigin);
+  }
+};
 
-  for (int y=0; y < coarseLattice.getNy(); ++y) {
-    T fEq[DESCRIPTOR<T>::q] {};
-    computeFeq(coarseLattice.get(x,y), fEq);
+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;
 
-    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      fineLattice.get(0,2*y)[iPop] = fEq[iPop] + coarse.getScalingFactor() * c2f_fneq[y][iPop];
-    }
+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:   " << 2*this->_coarseSize-1 << std::endl;
   }
 
-  for (int y=1; y < coarseLattice.getNy()-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];
+  void store()
+  {
+    auto& coarseLattice = this->_coarse.getSuperLattice();
 
-    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]
-                     );
+    for (int y=0; y < this->_coarseSize; ++y) {
+      const auto pos = this->getCoarseLatticeR(y);
 
-      fineLattice.get(0,1+2*y)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq;
+      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()
   {
-    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];
-
-    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(0,1)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq;
+    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()
   {
-    const int maxY = coarseLattice.getNy()-1;
-    const T rho = order3interpolation(
-                    c2f_rho[maxY  ],
-                    c2f_rho[maxY-1],
-                    c2f_rho[maxY-2]
-                  );
-    T u[2] {};
-    u[0] = order3interpolation(
-             c2f_u[maxY  ][0],
-             c2f_u[maxY-1][0],
-             c2f_u[maxY-2][0]
-           );
-    u[1] = order3interpolation(
-             c2f_u[maxY  ][1],
-             c2f_u[maxY-1][1],
-             c2f_u[maxY-2][1]
-           );
-    const T uSqr = u[0]*u[0] + u[1]*u[1];
-
-    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      const T fneq = order3interpolation(
-                       c2f_fneq[maxY  ][iPop],
-                       c2f_fneq[maxY-1][iPop],
-                       c2f_fneq[maxY-2][iPop]
-                     );
-      fineLattice.get(0,fineLattice.getNy()-2)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq;
+    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(2*(this->_coarseSize-1)-1);
+
+      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 BlockLatticeView2D<T,DESCRIPTOR>& lattice, int x, int y, T restrictedFneq[DESCRIPTOR<T>::q])
+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(x + DESCRIPTOR<T>::c[iPop][0], y + DESCRIPTOR<T>::c[iPop][1]), fNeq);
+    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];
@@ -513,24 +568,38 @@ void computeRestrictedFneq(const BlockLatticeView2D<T,DESCRIPTOR>& lattice, int
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
-void coupleF2C(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
-{
-  auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
-  auto& fineLattice   = fine.getSuperLattice().getBlockLattice(0);
+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();
 
-  const int x = coarseLattice.getNx()-1;
-  for (int y=0; y < coarseLattice.getNy(); ++y) {
-    T fEq[DESCRIPTOR<T>::q] {};
-    computeFeq(fineLattice.get(2,2*y), fEq);
+    for (int y=0; y < this->_coarseSize; ++y) {
+      const auto finePos   = this->getFineLatticeR(2*y);
+      const auto coarsePos = this->getCoarseLatticeR(y);
 
-    T fNeq[DESCRIPTOR<T>::q] {};
-    computeRestrictedFneq(fineLattice, 2, 2*y, fNeq);
+      T fEq[DESCRIPTOR<T>::q] {};
+      computeFeq(fineLattice.get(finePos), fEq);
 
-    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
-      coarseLattice.get(x,y)[iPop] = fEq[iPop] + coarse.getInvScalingFactor() * fNeq[iPop];
+      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[])
 {
@@ -543,6 +612,7 @@ int main(int argc, char* argv[])
   Vector<T,2> coarseOrigin { 0.0, 0.0 };
   Vector<T,2> coarseExtend { 0.5 * lx, ly  };
   IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
+
   Vector<T,2> fineOrigin   { 0.5 * lx - coarseDeltaX, 0.0 };
   Vector<T,2> fineExtend   { 0.5 * lx + coarseDeltaX, ly  };
   IndicatorCuboid2D<T> fineCuboid(fineExtend, fineOrigin);
@@ -592,10 +662,13 @@ int main(int argc, char* argv[])
     residuum);
   timer.start();
 
-  const auto sizeC2F = coarseGrid->getSuperLattice().getBlockLattice(0).getNy();
-  T c2f_rho[sizeC2F] {};
-  T c2f_u[sizeC2F][2] {};
-  T c2f_fneq[sizeC2F][DESCRIPTOR<T>::q] {};
+  const Vector<T,2> originC2F = fineOrigin;
+  const Vector<T,2> extendC2F { 0, fineExtend[1] };
+  FineCoupler<T,DESCRIPTOR> fineCoupler(*coarseGrid, *fineGrid, originC2F, extendC2F);
+
+  const Vector<T,2> originF2C { fineOrigin[0] + coarseDeltaX, fineOrigin[1] };
+  const Vector<T,2> extendF2C { 0, fineExtend[1] };
+  CoarseCoupler<T,DESCRIPTOR> coarseCoupler(*coarseGrid, *fineGrid, originF2C, extendF2C);
 
   for (int iT = 0; iT < coarseGrid->getConverter().getLatticeTime(maxPhysT); ++iT) {
     if (converge.hasConverged()) {
@@ -603,18 +676,18 @@ int main(int argc, char* argv[])
       break;
     }
 
-    storeC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq);
+    fineCoupler.store();
     coarseGrid->getSuperLattice().collideAndStream();
 
-    interpolateStoredC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq);
     fineGrid->getSuperLattice().collideAndStream();
-    coupleC2F(*coarseGrid, *fineGrid, c2f_rho, c2f_u, c2f_fneq);
+    fineCoupler.interpolate();
+    fineCoupler.couple();
 
-    storeC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq);
     fineGrid->getSuperLattice().collideAndStream();
-    coupleC2F(*coarseGrid, *fineGrid, c2f_rho, c2f_u, c2f_fneq);
+    fineCoupler.store();
+    fineCoupler.couple();
 
-    coupleF2C(*coarseGrid, *fineGrid);
+    coarseCoupler.couple();
 
     getResults(
       "coarse_",