summaryrefslogtreecommitdiff
path: root/apps/adrian/poiseuille2d
diff options
context:
space:
mode:
Diffstat (limited to 'apps/adrian/poiseuille2d')
-rw-r--r--apps/adrian/poiseuille2d/poiseuille2d.cpp80
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 };