diff options
| -rw-r--r-- | apps/adrian/poiseuille2d/poiseuille2d.cpp | 80 |
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 }; |