From 08e0d597bbf9c6f4c43df071410b69812ac9a8e2 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Tue, 15 Jan 2019 11:24:40 +0100 Subject: Interpolate vectors instead of scalars Same result, nicer code --- src/refinement/coupler2D.h | 4 +- src/refinement/coupler2D.hh | 129 ++++++++++++++------------------------------ 2 files changed, 41 insertions(+), 92 deletions(-) (limited to 'src/refinement') diff --git a/src/refinement/coupler2D.h b/src/refinement/coupler2D.h index db8ac15..d43a294 100644 --- a/src/refinement/coupler2D.h +++ b/src/refinement/coupler2D.h @@ -60,8 +60,8 @@ public: template class DESCRIPTOR> class FineCoupler2D : public Coupler2D { private: - std::vector _c2f_rho; - std::vector> _c2f_u; + std::vector> _c2f_rho; + std::vector::d>> _c2f_u; std::vector::q>> _c2f_fneq; public: 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::store() lbHelpers::computeRhoU(coarseCell, rho, u); lbHelpers::computeFneq(coarseCell, fNeq, rho, u); - _c2f_rho[y] = rho; + _c2f_rho[y] = Vector(rho); _c2f_u[y] = Vector(u); _c2f_fneq[y] = Vector::q>(fNeq); } } -template -T order4interpolation(T fm3, T fm1, T f1, T f3) +template +Vector order2interpolation(const Vector& f0, const Vector& f1) { - return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3); + return 0.5 * (f0 + f1); } -template -T order3interpolation(T fm1, T f1, T f3) +template +Vector order3interpolation(const std::vector>& 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 -T order2interpolation(T f0, T f1) +template +Vector order4interpolation(const std::vector>& 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 class DESCRIPTOR> void FineCoupler2D::interpolate() { @@ -153,15 +159,14 @@ void FineCoupler2D::interpolate() T rho{}; T u[2] {}; lbHelpers::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(rho), _c2f_rho[y]); + _c2f_u[y] = order2interpolation(Vector(u), _c2f_u[y]); T fNeq[DESCRIPTOR::q] {}; lbHelpers::computeFneq(coarseCell, fNeq, rho, u); - for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - _c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], _c2f_fneq[y][iPop]); - } + + _c2f_fneq[y] = order2interpolation(Vector::q>(fNeq), _c2f_fneq[y]); } } @@ -189,110 +194,54 @@ void FineCoupler2D::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 fineCell; fineLattice.get(finePos, fineCell); for (int iPop=0; iPop < DESCRIPTOR::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::equilibrium(iPop, rho, u, uSqr) + this->getScalingFactor() * fneq; + fineCell[iPop] = lbHelpers::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 fineCell; fineLattice.get(finePos, fineCell); for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - const T fneq = order3interpolation( - _c2f_fneq[0][iPop], - _c2f_fneq[1][iPop], - _c2f_fneq[2][iPop] - ); - fineCell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + this->getScalingFactor() * fneq; + fineCell[iPop] = lbHelpers::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 fineCell; fineLattice.get(finePos, fineCell); for (int iPop=0; iPop < DESCRIPTOR::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::equilibrium(iPop, rho, u, uSqr) + this->getScalingFactor() * fneq; + fineCell[iPop] = lbHelpers::equilibrium(iPop, rho[0], u.data, uSqr) + this->getScalingFactor() * fneq[iPop]; } fineLattice.set(finePos, fineCell); -- cgit v1.2.3