Basic grid refinement algorithm by Lagrava et al.

Starting point for integration into a more flexible setting.

See "Advances in multi-domain lattice Boltzmann grid refinement" (2012)

1 files changed, 189 insertions, 42 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 1808fbc..8bb982d 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -44,7 +44,7 @@ typedef double T;
 const T lx  = 2.;             // length of the channel
 const T ly  = 1.;             // height of the channel
 const int N = 16;             // resolution of the model
-const T Re = 10.;             // Reynolds number
+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
 const T residuum = 1e-5;      // residuum for the convergence check
@@ -274,6 +274,27 @@ public:
   }
 };
 
+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>
+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);
+}
+
 T order4interpolation(T fm3, T fm1, T f1, T f3)
 {
   return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3);
@@ -284,8 +305,73 @@ T order3interpolation(T fm1, T f1, T f3)
   return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3;
 }
 
+T order2interpolation(T f0, T f1)
+{
+  return 0.5 * (f0 + 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 i=0; i < DESCRIPTOR<T>::q; ++i) {
+      c2f_fneq[y][i] = fNeq[i];
+    }
+  }
+}
+
+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 i=0; i < DESCRIPTOR<T>::q; ++i) {
+      c2f_fneq[y][i] = order2interpolation(fNeq[i], c2f_fneq[y][i]);
+    }
+  }
+}
+
 template <typename T, template<typename> class DESCRIPTOR>
-void coupleC2F(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
+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);
@@ -293,54 +379,104 @@ void coupleC2F(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
   const int x = coarseLattice.getNx()-2;
 
   for (int y=0; y < coarseLattice.getNy(); ++y) {
+    T fEq[DESCRIPTOR<T>::q] {};
+    computeFeq(coarseLattice.get(x,y), fEq);
+
     for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      fineLattice.get(0, 2*y)[i] = coarseLattice.get(x,y)[i];
+      fineLattice.get(0,2*y)[i] = fEq[i] + coarse.getScalingFactor() * c2f_fneq[y][i];
     }
   }
 
   for (int y=2; 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];
+
     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]
-                                    );
+      const T fneqi = order4interpolation(
+                        c2f_fneq[y-1][i],
+                        c2f_fneq[y][i],
+                        c2f_fneq[y+1][i],
+                        c2f_fneq[y+2][i]
+                      );
+
+      fineLattice.get(0,1+2*y)[i] = lbHelpers<T,DESCRIPTOR>::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi;
     }
   }
 
-  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]
-                              );
-    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]
-        );
+  {
+    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 i=0; i < DESCRIPTOR<T>::q; ++i) {
+      const T fneqi = order3interpolation(
+                        c2f_fneq[0][i],
+                        c2f_fneq[1][i],
+                        c2f_fneq[2][i]
+                      );
+      fineLattice.get(0,1)[i] = lbHelpers<T,DESCRIPTOR>::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi;
+    }
   }
-}
 
-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);
-}
+  {
+    const T rho = order3interpolation(
+                    c2f_rho[coarseLattice.getNy()-1],
+                    c2f_rho[coarseLattice.getNy()-2],
+                    c2f_rho[coarseLattice.getNy()-3]
+                  );
+    T u[2] {};
+    u[0] = order3interpolation(
+             c2f_u[coarseLattice.getNy()-1][0],
+             c2f_u[coarseLattice.getNy()-2][0],
+             c2f_u[coarseLattice.getNy()-3][0]
+           );
+    u[1] = order3interpolation(
+             c2f_u[coarseLattice.getNy()-1][1],
+             c2f_u[coarseLattice.getNy()-2][1],
+             c2f_u[coarseLattice.getNy()-3][1]
+           );
+    const T uSqr = u[0]*u[0] + u[1]*u[1];
 
-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);
+    for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
+      const T fneqi = order3interpolation(
+                        c2f_fneq[coarseLattice.getNy()-1][i],
+                        c2f_fneq[coarseLattice.getNy()-2][i],
+                        c2f_fneq[coarseLattice.getNy()-3][i]
+                      );
+      fineLattice.get(0,fineLattice.getNy()-2)[i] = lbHelpers<T,DESCRIPTOR>::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi;
+    }
   }
 }
 
@@ -349,7 +485,7 @@ T computeRestrictedFneq(const BlockLatticeView2D<T,DESCRIPTOR>& lattice, int x,
 {
   T sum = T{0};
   for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-    T fNeq[DESCRIPTOR<T>::q]{};
+    T fNeq[DESCRIPTOR<T>::q] {};
     computeFneq(lattice.get(x + DESCRIPTOR<T>::c[i][0], y + DESCRIPTOR<T>::c[i][1]), fNeq);
     sum += fNeq[i];
   }
@@ -364,7 +500,7 @@ 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]{};
+    T fEq[DESCRIPTOR<T>::q] {};
     computeFeq(fineLattice.get(2,2*y), fEq);
 
     const T restrictedFneq = computeRestrictedFneq(fineLattice, 2, 2*y);
@@ -443,17 +579,28 @@ 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] {};
+
   for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
     if (converge.hasConverged()) {
       clout << "Simulation converged." << endl;
       break;
     }
 
+    storeC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq);
     coarseGrid.getSuperLattice().collideAndStream();
+    interpolateStoredC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq);
+
     fineGrid.getSuperLattice().collideAndStream();
-    coupleC2F(coarseGrid, fineGrid);
+    coupleC2F(coarseGrid, fineGrid, c2f_rho, c2f_u, c2f_fneq);
+
     fineGrid.getSuperLattice().collideAndStream();
-    coupleC2F(coarseGrid, fineGrid);
+    storeC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq);
+    coupleC2F(coarseGrid, fineGrid, c2f_rho, c2f_u, c2f_fneq);
+
     coupleF2C(coarseGrid, fineGrid);
 
     getResults(