diff options
Diffstat (limited to 'src/refinement')
| -rw-r--r-- | src/refinement/coupler2D.hh | 106 |
1 files changed, 43 insertions, 63 deletions
diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh index 49efe57..8370707 100644 --- a/src/refinement/coupler2D.hh +++ b/src/refinement/coupler2D.hh @@ -32,63 +32,6 @@ namespace olb { 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 iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { - fEq[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, 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); -} - -template <typename T> -T order4interpolation(T fm3, T fm1, T f1, T f3) -{ - return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3); -} - -template <typename T> -T order3interpolation(T fm1, T f1, T f3) -{ - return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3; -} - -template <typename T> -T order2interpolation(T f0, T f1) -{ - return 0.5 * (f0 + f1); -} - -template <typename T, template<typename> class DESCRIPTOR> -void computeRestrictedFneq(const SuperLattice2D<T,DESCRIPTOR>& lattice, Vector<int,3> pos, T restrictedFneq[DESCRIPTOR<T>::q]) -{ - for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { - T fNeq[DESCRIPTOR<T>::q] {}; - computeFneq(lattice.get(pos[0], pos[1] + DESCRIPTOR<T>::c[iPop][0], pos[2] + DESCRIPTOR<T>::c[iPop][1]), fNeq); - - for (int jPop=0; jPop < DESCRIPTOR<T>::q; ++jPop) { - restrictedFneq[jPop] += fNeq[jPop]; - } - } - - for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { - restrictedFneq[iPop] /= DESCRIPTOR<T>::q; - } -} - - -template <typename T, template<typename> class DESCRIPTOR> Vector<int,3> Coupler2D<T,DESCRIPTOR>::getFineLatticeR(int y) const { if (_vertical) { @@ -159,7 +102,7 @@ void FineCoupler2D<T,DESCRIPTOR>::store() T u[2] {}; T fNeq[DESCRIPTOR<T>::q] {}; coarseLattice.get(pos).computeRhoU(rho, u); - computeFneq(coarseLattice.get(pos), fNeq); + coarseLattice.get(pos).computeFneq(fNeq); _c2f_rho[y] = rho; _c2f_u[y] = Vector<T,2>(u); @@ -167,23 +110,41 @@ void FineCoupler2D<T,DESCRIPTOR>::store() } } +template <typename T> +T order4interpolation(T fm3, T fm1, T f1, T f3) +{ + return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3); +} + +template <typename T> +T order3interpolation(T fm1, T f1, T f3) +{ + return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3; +} + +template <typename T> +T order2interpolation(T f0, T f1) +{ + return 0.5 * (f0 + f1); +} + template <typename T, template<typename> class DESCRIPTOR> void FineCoupler2D<T,DESCRIPTOR>::interpolate() { auto& coarseLattice = this->_coarse.getSuperLattice(); for (int y=0; y < this->_coarseSize; ++y) { - const auto coarsePos = this->getCoarseLatticeR(y); + const auto coarseCell = coarseLattice.get(this->getCoarseLatticeR(y)); T rho{}; T u[2] {}; - coarseLattice.get(coarsePos).computeRhoU(rho, u); + coarseCell.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(coarsePos), fNeq); + coarseCell.computeFneq(fNeq); for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { _c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], _c2f_fneq[y][iPop]); } @@ -201,7 +162,7 @@ void FineCoupler2D<T,DESCRIPTOR>::couple() const auto finePos = this->getFineLatticeR(2*y); T fEq[DESCRIPTOR<T>::q] {}; - computeFeq(coarseLattice.get(coarsePos), fEq); + coarseLattice.get(coarsePos).computeFeq(fEq); for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { fineLattice.get(finePos)[iPop] = fEq[iPop] + this->_coarse.getScalingFactor() * _c2f_fneq[y][iPop]; @@ -309,6 +270,25 @@ void FineCoupler2D<T,DESCRIPTOR>::couple() template <typename T, template<typename> class DESCRIPTOR> +void computeRestrictedFneq(const SuperLattice2D<T,DESCRIPTOR>& lattice, + Vector<int,3> pos, + T restrictedFneq[DESCRIPTOR<T>::q]) +{ + for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { + T fNeq[DESCRIPTOR<T>::q] {}; + lattice.get(pos[0], pos[1] + DESCRIPTOR<T>::c[iPop][0], pos[2] + DESCRIPTOR<T>::c[iPop][1]).computeFneq(fNeq); + + for (int jPop=0; jPop < DESCRIPTOR<T>::q; ++jPop) { + restrictedFneq[jPop] += fNeq[jPop]; + } + } + + for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { + restrictedFneq[iPop] /= DESCRIPTOR<T>::q; + } +} + +template <typename T, template<typename> class DESCRIPTOR> CoarseCoupler2D<T,DESCRIPTOR>::CoarseCoupler2D( Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRIPTOR>& fine, Vector<T,2> origin, Vector<T,2> extend): @@ -331,7 +311,7 @@ void CoarseCoupler2D<T,DESCRIPTOR>::couple() const auto coarsePos = this->getCoarseLatticeR(y); T fEq[DESCRIPTOR<T>::q] {}; - computeFeq(fineLattice.get(finePos), fEq); + fineLattice.get(finePos).computeFeq(fEq); T fNeq[DESCRIPTOR<T>::q] {}; computeRestrictedFneq(fineLattice, finePos, fNeq); |