1 files changed, 64 insertions, 16 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 583282e..1808fbc 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -43,7 +43,7 @@ typedef double T;
 
 const T lx  = 2.;             // length of the channel
 const T ly  = 1.;             // height of the channel
-const int N = 20;             // resolution of the model
+const int N = 16;             // resolution of the model
 const T Re = 10.;             // Reynolds number
 const T maxPhysT = 40.;       // max. simulation time in s, SI unit
 const T physInterval = 0.25;  // interval for the convergence check in s
@@ -262,6 +262,16 @@ public:
   {
     return *_lattice;
   }
+
+  T getScalingFactor()
+  {
+    return 4. - _converter->getLatticeRelaxationFrequency();
+  }
+
+  T getInvScalingFactor()
+  {
+    return 1./getScalingFactor();
+  }
 };
 
 T order4interpolation(T fm3, T fm1, T f1, T f3)
@@ -291,26 +301,59 @@ void coupleC2F(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
   for (int y=2; y < coarseLattice.getNy()-2; ++y) {
     for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
       fineLattice.get(0,1+2*y)[i] = order4interpolation(
-        coarseLattice.get(x,y-1)[i],
-        coarseLattice.get(x,y)[i],
-        coarseLattice.get(x,y+1)[i],
-        coarseLattice.get(x,y+2)[i]
-      );
+                                      coarseLattice.get(x,y-1)[i],
+                                      coarseLattice.get(x,y)[i],
+                                      coarseLattice.get(x,y+1)[i],
+                                      coarseLattice.get(x,y+2)[i]
+                                    );
     }
   }
 
   for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
     fineLattice.get(0,1)[i] = order3interpolation(
-      coarseLattice.get(x,0)[i],
-      coarseLattice.get(x,1)[i],
-      coarseLattice.get(x,2)[i]
-    );
+                                coarseLattice.get(x,0)[i],
+                                coarseLattice.get(x,1)[i],
+                                coarseLattice.get(x,2)[i]
+                              );
     fineLattice.get(0,fineLattice.getNy()-2)[i] = order3interpolation(
-      coarseLattice.get(x,coarseLattice.getNy()-1)[i],
-      coarseLattice.get(x,coarseLattice.getNy()-2)[i],
-      coarseLattice.get(x,coarseLattice.getNy()-3)[i]
-    );
+          coarseLattice.get(x,coarseLattice.getNy()-1)[i],
+          coarseLattice.get(x,coarseLattice.getNy()-2)[i],
+          coarseLattice.get(x,coarseLattice.getNy()-3)[i]
+        );
+  }
+}
+
+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, 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 i=0; i < DESCRIPTOR<T>::q; ++i) {
+    fEq[i] = lbHelpers<T,DESCRIPTOR>::equilibrium(i, rho, u, uSqr);
+  }
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+T computeRestrictedFneq(const BlockLatticeView2D<T,DESCRIPTOR>& lattice, int x, int y)
+{
+  T sum = T{0};
+  for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
+    T fNeq[DESCRIPTOR<T>::q]{};
+    computeFneq(lattice.get(x + DESCRIPTOR<T>::c[i][0], y + DESCRIPTOR<T>::c[i][1]), fNeq);
+    sum += fNeq[i];
   }
+  return sum / DESCRIPTOR<T>::q;
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
@@ -321,8 +364,13 @@ void coupleF2C(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
 
   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);
+
+    const T restrictedFneq = computeRestrictedFneq(fineLattice, 2, 2*y);
+
     for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      coarseLattice.get(x,y)[i] = fineLattice.get(2,2*y)[i];
+      coarseLattice.get(x,y)[i] = fEq[i] + coarse.getInvScalingFactor() * restrictedFneq;
     }
   }
 }
@@ -335,7 +383,7 @@ int main(int argc, char* argv[])
 
   const T coarseDeltaX = 1./N;
 
-  Vector<T,2> coarseOrigin { 0.0     , 0.0 };
+  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 };