diff options
Diffstat (limited to 'src/refinement/coupler2D.hh')
| -rw-r--r-- | src/refinement/coupler2D.hh | 129 |
1 files changed, 39 insertions, 90 deletions
diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh index 1726df4..2218c28 100644 --- a/src/refinement/coupler2D.hh +++ b/src/refinement/coupler2D.hh @@ -117,30 +117,36 @@ void FineCoupler2D<T,DESCRIPTOR>::store() lbHelpers<T,DESCRIPTOR>::computeRhoU(coarseCell, rho, u); lbHelpers<T,DESCRIPTOR>::computeFneq(coarseCell, fNeq, rho, u); - _c2f_rho[y] = rho; + _c2f_rho[y] = Vector<T,1>(rho); _c2f_u[y] = Vector<T,2>(u); _c2f_fneq[y] = Vector<T,DESCRIPTOR<T>::q>(fNeq); } } -template <typename T> -T order4interpolation(T fm3, T fm1, T f1, T f3) +template <typename T, unsigned N> +Vector<T,N> order2interpolation(const Vector<T,N>& f0, const Vector<T,N>& f1) { - return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3); + return 0.5 * (f0 + f1); } -template <typename T> -T order3interpolation(T fm1, T f1, T f3) +template <typename T, unsigned N> +Vector<T,N> order3interpolation(const std::vector<Vector<T,N>>& data, int y, bool ascending) { - return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3; + if (ascending) { + return 3./8. * data[y] + 3./4. * data[y+1] - 1./8. * data[y+2]; + } + else { + return 3./8. * data[y] + 3./4. * data[y-1] - 1./8. * data[y-2]; + } } -template <typename T> -T order2interpolation(T f0, T f1) +template <typename T, unsigned N> +Vector<T,N> order4interpolation(const std::vector<Vector<T,N>>& data, int y) { - return 0.5 * (f0 + f1); + return 9./16. * (data[y] + data[y+1]) - 1./16. * (data[y-1] + data[y+2]); } + template <typename T, template<typename> class DESCRIPTOR> void FineCoupler2D<T,DESCRIPTOR>::interpolate() { @@ -153,15 +159,14 @@ void FineCoupler2D<T,DESCRIPTOR>::interpolate() T rho{}; T u[2] {}; lbHelpers<T,DESCRIPTOR>::computeRhoU(coarseCell, 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]); + + _c2f_rho[y] = order2interpolation(Vector<T,1>(rho), _c2f_rho[y]); + _c2f_u[y] = order2interpolation(Vector<T,2>(u), _c2f_u[y]); T fNeq[DESCRIPTOR<T>::q] {}; lbHelpers<T,DESCRIPTOR>::computeFneq(coarseCell, fNeq, rho, u); - for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { - _c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], _c2f_fneq[y][iPop]); - } + + _c2f_fneq[y] = order2interpolation(Vector<T,DESCRIPTOR<T>::q>(fNeq), _c2f_fneq[y]); } } @@ -189,110 +194,54 @@ void FineCoupler2D<T,DESCRIPTOR>::couple() } for (int y=1; y < this->_coarseSize-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]; + const auto rho = order4interpolation(_c2f_rho, y); + const auto u = order4interpolation(_c2f_u, y); + const auto fneq = order4interpolation(_c2f_fneq, y); + + const T uSqr = u*u; const auto finePos = this->getFineLatticeR(1+2*y); Cell<T,DESCRIPTOR> fineCell; fineLattice.get(finePos, fineCell); for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { - const T fneq = order4interpolation( - _c2f_fneq[y-1][iPop], - _c2f_fneq[y][iPop], - _c2f_fneq[y+1][iPop], - _c2f_fneq[y+2][iPop] - ); - - fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->getScalingFactor() * fneq; + fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop]; } fineLattice.set(finePos, fineCell); } { - 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]; + const auto rho = order3interpolation(_c2f_rho, 0, true); + const auto u = order3interpolation(_c2f_u, 0, true); + const auto fneq = order3interpolation(_c2f_fneq, 0, true); + + const T uSqr = u*u; const auto& finePos = this->getFineLatticeR(1); Cell<T,DESCRIPTOR> fineCell; fineLattice.get(finePos, fineCell); for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { - const T fneq = order3interpolation( - _c2f_fneq[0][iPop], - _c2f_fneq[1][iPop], - _c2f_fneq[2][iPop] - ); - fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->getScalingFactor() * fneq; + fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop]; } fineLattice.set(finePos, fineCell); } { - const T rho = order3interpolation( - _c2f_rho[this->_coarseSize-1], - _c2f_rho[this->_coarseSize-2], - _c2f_rho[this->_coarseSize-3] - ); - T u[2] {}; - u[0] = order3interpolation( - _c2f_u[this->_coarseSize-1][0], - _c2f_u[this->_coarseSize-2][0], - _c2f_u[this->_coarseSize-3][0] - ); - u[1] = order3interpolation( - _c2f_u[this->_coarseSize-1][1], - _c2f_u[this->_coarseSize-2][1], - _c2f_u[this->_coarseSize-3][1] - ); - const T uSqr = u[0]*u[0] + u[1]*u[1]; + const auto rho = order3interpolation(_c2f_rho, this->_coarseSize-1, false); + const auto u = order3interpolation(_c2f_u, this->_coarseSize-1, false); + const auto fneq = order3interpolation(_c2f_fneq, this->_coarseSize-1, false); + + const T uSqr = u*u; const auto& finePos = this->getFineLatticeR(this->_fineSize-2); Cell<T,DESCRIPTOR> fineCell; fineLattice.get(finePos, fineCell); for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) { - const T fneq = order3interpolation( - _c2f_fneq[this->_coarseSize-1][iPop], - _c2f_fneq[this->_coarseSize-2][iPop], - _c2f_fneq[this->_coarseSize-3][iPop] - ); - fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->getScalingFactor() * fneq; + fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop]; } fineLattice.set(finePos, fineCell); |