From 663d1dbbdc47e6d2863eb7b397146b3f41795073 Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Wed, 2 Jan 2019 12:41:52 +0100
Subject: 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)
---
apps/adrian/poiseuille2d/poiseuille2d.cpp | 231 ++++++++++++++++++++++++------
1 file changed, 189 insertions(+), 42 deletions(-)
(limited to 'apps/adrian')
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(
--
cgit v1.2.3