summaryrefslogtreecommitdiff
path: root/src/refinement
diff options
context:
space:
mode:
Diffstat (limited to 'src/refinement')
-rw-r--r--src/refinement/coupler2D.h4
-rw-r--r--src/refinement/coupler2D.hh129
2 files changed, 41 insertions, 92 deletions
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 <typename T, template<typename> class DESCRIPTOR>
class FineCoupler2D : public Coupler2D<T,DESCRIPTOR> {
private:
- std::vector<T> _c2f_rho;
- std::vector<Vector<T,2>> _c2f_u;
+ std::vector<Vector<T,1>> _c2f_rho;
+ std::vector<Vector<T,DESCRIPTOR<T>::d>> _c2f_u;
std::vector<Vector<T,DESCRIPTOR<T>::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<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);