Move distribution decomposition methods to cell

1 files changed, 43 insertions, 63 deletions

diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh
index 49efe57..8370707 100644
--- a/src/refinement/coupler2D.hh
+++ b/src/refinement/coupler2D.hh
@@ -32,63 +32,6 @@ namespace olb {
 
 
 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);
-}
-
-template <typename T>
-T order4interpolation(T fm3, T fm1, T f1, T f3)
-{
-  return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3);
-}
-
-template <typename T>
-T order3interpolation(T fm1, T f1, T f3)
-{
-  return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3;
-}
-
-template <typename T>
-T order2interpolation(T f0, T f1)
-{
-  return 0.5 * (f0 + f1);
-}
-
-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>
 Vector<int,3> Coupler2D<T,DESCRIPTOR>::getFineLatticeR(int y) const
 {
   if (_vertical) {
@@ -159,7 +102,7 @@ void FineCoupler2D<T,DESCRIPTOR>::store()
     T u[2] {};
     T fNeq[DESCRIPTOR<T>::q] {};
     coarseLattice.get(pos).computeRhoU(rho, u);
-    computeFneq(coarseLattice.get(pos), fNeq);
+    coarseLattice.get(pos).computeFneq(fNeq);
 
     _c2f_rho[y]  = rho;
     _c2f_u[y]    = Vector<T,2>(u);
@@ -167,23 +110,41 @@ void FineCoupler2D<T,DESCRIPTOR>::store()
   }
 }
 
+template <typename T>
+T order4interpolation(T fm3, T fm1, T f1, T f3)
+{
+  return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3);
+}
+
+template <typename T>
+T order3interpolation(T fm1, T f1, T f3)
+{
+  return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3;
+}
+
+template <typename T>
+T order2interpolation(T f0, T f1)
+{
+  return 0.5 * (f0 + f1);
+}
+
 template <typename T, template<typename> class DESCRIPTOR>
 void FineCoupler2D<T,DESCRIPTOR>::interpolate()
 {
   auto& coarseLattice = this->_coarse.getSuperLattice();
 
   for (int y=0; y < this->_coarseSize; ++y) {
-    const auto coarsePos = this->getCoarseLatticeR(y);
+    const auto coarseCell = coarseLattice.get(this->getCoarseLatticeR(y));
 
     T rho{};
     T u[2] {};
-    coarseLattice.get(coarsePos).computeRhoU(rho, u);
+    coarseCell.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);
+    coarseCell.computeFneq(fNeq);
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
       _c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], _c2f_fneq[y][iPop]);
     }
@@ -201,7 +162,7 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
     const auto finePos   = this->getFineLatticeR(2*y);
 
     T fEq[DESCRIPTOR<T>::q] {};
-    computeFeq(coarseLattice.get(coarsePos), fEq);
+    coarseLattice.get(coarsePos).computeFeq(fEq);
 
     for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
       fineLattice.get(finePos)[iPop] = fEq[iPop] + this->_coarse.getScalingFactor() * _c2f_fneq[y][iPop];
@@ -309,6 +270,25 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
 
 
 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] {};
+    lattice.get(pos[0], pos[1] + DESCRIPTOR<T>::c[iPop][0], pos[2] + DESCRIPTOR<T>::c[iPop][1]).computeFneq(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>
 CoarseCoupler2D<T,DESCRIPTOR>::CoarseCoupler2D(
   Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRIPTOR>& fine,
   Vector<T,2> origin, Vector<T,2> extend):
@@ -331,7 +311,7 @@ void CoarseCoupler2D<T,DESCRIPTOR>::couple()
     const auto coarsePos = this->getCoarseLatticeR(y);
 
     T fEq[DESCRIPTOR<T>::q] {};
-    computeFeq(fineLattice.get(finePos), fEq);
+    fineLattice.get(finePos).computeFeq(fEq);
 
     T fNeq[DESCRIPTOR<T>::q] {};
     computeRestrictedFneq(fineLattice, finePos, fNeq);