Verify mass loss when using linear interpolation for C2F

2 files changed, 13 insertions, 7 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 19dc4d7..b49f24c 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -36,10 +36,10 @@ typedef double T;
 
 #define DESCRIPTOR descriptors::D2Q9Descriptor
 
-const T lx  = 4.0;            // length of the channel
-const T ly  = 1.0;            // height of the channel
-const int N = 50;             // resolution of the model
-const T Re  = 100.;           // Reynolds number
+const T lx   = 4.0;            // length of the channel
+const T ly   = 1.0;            // height of the channel
+const int N  = 50;             // resolution of the model
+const T Re   = 100.;           // Reynolds number
 const T uMax = 0.02;          // Max lattice speed
 const T maxPhysT = 60.;       // max. simulation time in s, SI unit
 const T physInterval = 0.25;  // interval for the convergence check in s
diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh
index 91a7db8..78aad46 100644
--- a/src/refinement/coupler2D.hh
+++ b/src/refinement/coupler2D.hh
@@ -131,6 +131,12 @@ Vector<T,N> order2interpolation(const Vector<T,N>& f0, const Vector<T,N>& f1)
 }
 
 template <typename T, unsigned N>
+Vector<T,N> order2interpolation(const std::vector<Vector<T,N>>& data, int y)
+{
+  return 0.5 * (data[y] + data[y+1]);
+}
+
+template <typename T, unsigned N>
 Vector<T,N> order3interpolation(const std::vector<Vector<T,N>>& data, int y, bool ascending)
 {
   if (ascending) {
@@ -195,9 +201,9 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
   }
 
   for (int y=1; y < this->_coarseSize-2; ++y) {
-    const auto rho  = order4interpolation(_c2f_rho,  y);
-    const auto u    = order4interpolation(_c2f_u,    y);
-    const auto fneq = order4interpolation(_c2f_fneq, y);
+    const auto rho  = order2interpolation(_c2f_rho,  y);
+    const auto u    = order2interpolation(_c2f_u,    y);
+    const auto fneq = order2interpolation(_c2f_fneq, y);
 
     const T uSqr = u*u;