summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/refinement/coupler2D.hh106
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);