summaryrefslogtreecommitdiff
path: root/apps/adrian
diff options
context:
space:
mode:
authorAdrian Kummerlaender2019-01-05 17:45:46 +0100
committerAdrian Kummerlaender2019-06-24 15:14:37 +0200
commit6b501a6981f001cb0b20404b51fe9ec6eb6f4325 (patch)
treed5e04532929899a3bd43c40ae7b666a73a5ff95c /apps/adrian
parentd0f1e440eaa92116781f69c953db47f758bde9b4 (diff)
downloadgrid_refinement_openlb-6b501a6981f001cb0b20404b51fe9ec6eb6f4325.tar
grid_refinement_openlb-6b501a6981f001cb0b20404b51fe9ec6eb6f4325.tar.gz
grid_refinement_openlb-6b501a6981f001cb0b20404b51fe9ec6eb6f4325.tar.bz2
grid_refinement_openlb-6b501a6981f001cb0b20404b51fe9ec6eb6f4325.tar.lz
grid_refinement_openlb-6b501a6981f001cb0b20404b51fe9ec6eb6f4325.tar.xz
grid_refinement_openlb-6b501a6981f001cb0b20404b51fe9ec6eb6f4325.tar.zst
grid_refinement_openlb-6b501a6981f001cb0b20404b51fe9ec6eb6f4325.zip
Refactor grid coupling into classes
Diffstat (limited to 'apps/adrian')
-rw-r--r--apps/adrian/poiseuille2d/poiseuille2d.cpp417
1 files changed, 245 insertions, 172 deletions
diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 22c35cb..632fba5 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -38,7 +38,7 @@ typedef double T;
const T lx = 4.; // length of the channel
const T ly = 1.; // height of the channel
const int N = 30; // resolution of the model
-const T Re = 100.; // Reynolds number
+const T Re = 100.; // Reynolds number
const T maxPhysT = 60.; // max. simulation time in s, SI unit
const T physInterval = 0.25; // interval for the convergence check in s
const T residuum = 1e-5; // residuum for the convergence check
@@ -325,182 +325,237 @@ T order2interpolation(T f0, T f1)
}
template <typename T, template<typename> class DESCRIPTOR>
-void storeC2F(
- Grid<T,DESCRIPTOR>& coarse,
- T c2f_rho[],
- T c2f_u[][2],
- T c2f_fneq[][DESCRIPTOR<T>::q]
-)
-{
- auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
-
- const int x = coarseLattice.getNx()-2;
- for (int y=0; y < coarseLattice.getNy(); ++y) {
- T rho{};
- T u[2] {};
- coarseLattice.get(x,y).computeRhoU(rho, u);
- c2f_rho[y] = rho;
- c2f_u[y][0] = u[0];
- c2f_u[y][1] = u[1];
-
- T fNeq[DESCRIPTOR<T>::q] {};
- computeFneq(coarseLattice.get(x,y), fNeq);
- for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
- c2f_fneq[y][iPop] = fNeq[iPop];
+class Coupler {
+protected:
+ Grid<T,DESCRIPTOR>& _coarse;
+ Grid<T,DESCRIPTOR>& _fine;
+ const Vector<T,2> _origin;
+ const Vector<T,2> _extend;
+ const int _coarseSize;
+ const bool _vertical;
+ Vector<int,3> _coarseOrigin;
+ Vector<int,3> _fineOrigin;
+
+ Vector<int,3> getFineLatticeR(int y) const
+ {
+ if (_vertical) {
+ return _fineOrigin + Vector<int,3> {0, 0, y};
+ }
+ else {
+ return _fineOrigin + Vector<int,3> {0, y, 0};
}
}
-}
-template <typename T, template<typename> class DESCRIPTOR>
-void interpolateStoredC2F(
- Grid<T,DESCRIPTOR>& coarse,
- T c2f_rho[],
- T c2f_u[][2],
- T c2f_fneq[][DESCRIPTOR<T>::q]
-)
-{
- auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
-
- const int x = coarseLattice.getNx()-2;
- for (int y=0; y < coarseLattice.getNy(); ++y) {
- T rho{};
- T u[2] {};
- coarseLattice.get(x,y).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(x,y), fNeq);
- for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
- c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], c2f_fneq[y][iPop]);
+ Vector<int,3> getCoarseLatticeR(int y) const
+ {
+ if (_vertical) {
+ return _coarseOrigin + Vector<int,3> {0, 0, y};
+ }
+ else {
+ return _coarseOrigin + Vector<int,3> {0, y, 0};
}
}
-}
-template <typename T, template<typename> class DESCRIPTOR>
-void coupleC2F(
- Grid<T,DESCRIPTOR>& coarse,
- Grid<T,DESCRIPTOR>& fine,
- const T c2f_rho[],
- const T c2f_u[][2],
- const T c2f_fneq[][DESCRIPTOR<T>::q]
-)
-{
- auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
- auto& fineLattice = fine.getSuperLattice().getBlockLattice(0);
+public:
+ Coupler(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine, Vector<T,2> origin, Vector<T,2> extend):
+ _coarse(coarse),
+ _fine(fine),
+ _origin(origin),
+ _extend(extend),
+ _coarseSize(extend.norm() / coarse.getConverter().getPhysDeltaX()),
+ _vertical(util::nearZero(extend[0]))
+ {
+ OLB_ASSERT(util::nearZero(extend[0]) || util::nearZero(extend[1]), "Coupling is only implemented alongside unit vectors");
- const int x = coarseLattice.getNx()-2;
+ _coarse.getCuboidGeometry().getFloorLatticeR(_origin, _coarseOrigin);
+ _fine.getCuboidGeometry().getFloorLatticeR(_origin, _fineOrigin);
+ }
+};
- for (int y=0; y < coarseLattice.getNy(); ++y) {
- T fEq[DESCRIPTOR<T>::q] {};
- computeFeq(coarseLattice.get(x,y), fEq);
+template <typename T, template<typename> class DESCRIPTOR>
+class FineCoupler : public Coupler<T,DESCRIPTOR> {
+private:
+ std::vector<T> _c2f_rho;
+ std::vector<Vector<T,2>> _c2f_u;
+ std::vector<Vector<T,DESCRIPTOR<T>::q>> _c2f_fneq;
- for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
- fineLattice.get(0,2*y)[iPop] = fEq[iPop] + coarse.getScalingFactor() * c2f_fneq[y][iPop];
- }
+public:
+ FineCoupler(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine, Vector<T,2> origin, Vector<T,2> extend):
+ Coupler<T,DESCRIPTOR>(coarse, fine, origin, extend),
+ _c2f_rho(this->_coarseSize),
+ _c2f_u(this->_coarseSize, Vector<T,2>(T{})),
+ _c2f_fneq(this->_coarseSize, Vector<T,DESCRIPTOR<T>::q>(T{}))
+ {
+ OstreamManager clout(std::cout,"C2F");
+ clout << "coarse origin: " << this->_coarseOrigin[0] << " " << this->_coarseOrigin[1] << " " << this->_coarseOrigin[2] << std::endl;
+ clout << "fine origin: " << this->_fineOrigin[0] << " " << this->_fineOrigin[1] << " " << this->_fineOrigin[2] << std::endl;
+ clout << "fine size: " << 2*this->_coarseSize-1 << std::endl;
}
- for (int y=1; y < coarseLattice.getNy()-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];
+ void store()
+ {
+ auto& coarseLattice = this->_coarse.getSuperLattice();
- 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]
- );
+ for (int y=0; y < this->_coarseSize; ++y) {
+ const auto pos = this->getCoarseLatticeR(y);
- fineLattice.get(0,1+2*y)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq;
+ T rho{};
+ T u[2] {};
+ T fNeq[DESCRIPTOR<T>::q] {};
+ coarseLattice.get(pos).computeRhoU(rho, u);
+ computeFneq(coarseLattice.get(pos), fNeq);
+
+ _c2f_rho[y] = rho;
+ _c2f_u[y] = Vector<T,2>(u);
+ _c2f_fneq[y] = Vector<T,DESCRIPTOR<T>::q>(fNeq);
}
}
+ void interpolate()
{
- 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];
-
- 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]
- );
- fineLattice.get(0,1)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq;
+ auto& coarseLattice = this->_coarse.getSuperLattice();
+
+ for (int y=0; y < this->_coarseSize; ++y) {
+ const auto coarsePos = this->getCoarseLatticeR(y);
+
+ T rho{};
+ T u[2] {};
+ coarseLattice.get(coarsePos).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);
+ for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+ _c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], _c2f_fneq[y][iPop]);
+ }
}
}
+ void couple()
{
- const int maxY = coarseLattice.getNy()-1;
- const T rho = order3interpolation(
- c2f_rho[maxY ],
- c2f_rho[maxY-1],
- c2f_rho[maxY-2]
- );
- T u[2] {};
- u[0] = order3interpolation(
- c2f_u[maxY ][0],
- c2f_u[maxY-1][0],
- c2f_u[maxY-2][0]
- );
- u[1] = order3interpolation(
- c2f_u[maxY ][1],
- c2f_u[maxY-1][1],
- c2f_u[maxY-2][1]
- );
- const T uSqr = u[0]*u[0] + u[1]*u[1];
-
- for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
- const T fneq = order3interpolation(
- c2f_fneq[maxY ][iPop],
- c2f_fneq[maxY-1][iPop],
- c2f_fneq[maxY-2][iPop]
- );
- fineLattice.get(0,fineLattice.getNy()-2)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq;
+ auto& coarseLattice = this->_coarse.getSuperLattice();
+ auto& fineLattice = this->_fine.getSuperLattice();
+
+ for (int y=0; y < this->_coarseSize; ++y) {
+ const auto coarsePos = this->getCoarseLatticeR(y);
+ const auto finePos = this->getFineLatticeR(2*y);
+
+ T fEq[DESCRIPTOR<T>::q] {};
+ computeFeq(coarseLattice.get(coarsePos), fEq);
+
+ for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+ fineLattice.get(finePos)[iPop] = fEq[iPop] + this->_coarse.getScalingFactor() * _c2f_fneq[y][iPop];
+ }
+ }
+
+ 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 finePos = this->getFineLatticeR(1+2*y);
+
+ 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]
+ );
+
+ fineLattice.get(finePos)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
+ }
+ }
+
+ {
+ 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 finePos = this->getFineLatticeR(1);
+
+ 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]
+ );
+ fineLattice.get(finePos)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
+ }
+ }
+
+ {
+ 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 finePos = this->getFineLatticeR(2*(this->_coarseSize-1)-1);
+
+ 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]
+ );
+ fineLattice.get(finePos)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq;
+ }
}
}
-}
+};
template <typename T, template<typename> class DESCRIPTOR>
-void computeRestrictedFneq(const BlockLatticeView2D<T,DESCRIPTOR>& lattice, int x, int y, T restrictedFneq[DESCRIPTOR<T>::q])
+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(x + DESCRIPTOR<T>::c[iPop][0], y + DESCRIPTOR<T>::c[iPop][1]), fNeq);
+ 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];
@@ -513,24 +568,38 @@ void computeRestrictedFneq(const BlockLatticeView2D<T,DESCRIPTOR>& lattice, int
}
template <typename T, template<typename> class DESCRIPTOR>
-void coupleF2C(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
-{
- auto& coarseLattice = coarse.getSuperLattice().getBlockLattice(0);
- auto& fineLattice = fine.getSuperLattice().getBlockLattice(0);
+class CoarseCoupler : public Coupler<T,DESCRIPTOR> {
+public:
+ CoarseCoupler(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine, Vector<T,2> origin, Vector<T,2> extend):
+ Coupler<T,DESCRIPTOR>(coarse, fine, origin, extend)
+ {
+ OstreamManager clout(std::cout,"F2C");
+ clout << "coarse origin: " << this->_coarseOrigin[0] << " " << this->_coarseOrigin[1] << " " << this->_coarseOrigin[2] << std::endl;
+ clout << "fine origin: " << this->_fineOrigin[0] << " " << this->_fineOrigin[1] << " " << this->_fineOrigin[2] << std::endl;
+ clout << "coarse size: " << this->_coarseSize << std::endl;
+ }
+
+ void couple()
+ {
+ auto& fineLattice = this->_fine.getSuperLattice();
+ auto& coarseLattice = this->_coarse.getSuperLattice();
- const int x = coarseLattice.getNx()-1;
- for (int y=0; y < coarseLattice.getNy(); ++y) {
- T fEq[DESCRIPTOR<T>::q] {};
- computeFeq(fineLattice.get(2,2*y), fEq);
+ for (int y=0; y < this->_coarseSize; ++y) {
+ const auto finePos = this->getFineLatticeR(2*y);
+ const auto coarsePos = this->getCoarseLatticeR(y);
- T fNeq[DESCRIPTOR<T>::q] {};
- computeRestrictedFneq(fineLattice, 2, 2*y, fNeq);
+ T fEq[DESCRIPTOR<T>::q] {};
+ computeFeq(fineLattice.get(finePos), fEq);
- for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
- coarseLattice.get(x,y)[iPop] = fEq[iPop] + coarse.getInvScalingFactor() * fNeq[iPop];
+ T fNeq[DESCRIPTOR<T>::q] {};
+ computeRestrictedFneq(fineLattice, finePos, fNeq);
+
+ for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+ coarseLattice.get(coarsePos)[iPop] = fEq[iPop] + this->_coarse.getInvScalingFactor() * fNeq[iPop];
+ }
}
}
-}
+};
int main(int argc, char* argv[])
{
@@ -543,6 +612,7 @@ int main(int argc, char* argv[])
Vector<T,2> coarseOrigin { 0.0, 0.0 };
Vector<T,2> coarseExtend { 0.5 * lx, ly };
IndicatorCuboid2D<T> coarseCuboid(coarseExtend, coarseOrigin);
+
Vector<T,2> fineOrigin { 0.5 * lx - coarseDeltaX, 0.0 };
Vector<T,2> fineExtend { 0.5 * lx + coarseDeltaX, ly };
IndicatorCuboid2D<T> fineCuboid(fineExtend, fineOrigin);
@@ -592,10 +662,13 @@ int main(int argc, char* argv[])
residuum);
timer.start();
- const auto sizeC2F = coarseGrid->getSuperLattice().getBlockLattice(0).getNy();
- T c2f_rho[sizeC2F] {};
- T c2f_u[sizeC2F][2] {};
- T c2f_fneq[sizeC2F][DESCRIPTOR<T>::q] {};
+ const Vector<T,2> originC2F = fineOrigin;
+ const Vector<T,2> extendC2F { 0, fineExtend[1] };
+ FineCoupler<T,DESCRIPTOR> fineCoupler(*coarseGrid, *fineGrid, originC2F, extendC2F);
+
+ const Vector<T,2> originF2C { fineOrigin[0] + coarseDeltaX, fineOrigin[1] };
+ const Vector<T,2> extendF2C { 0, fineExtend[1] };
+ CoarseCoupler<T,DESCRIPTOR> coarseCoupler(*coarseGrid, *fineGrid, originF2C, extendF2C);
for (int iT = 0; iT < coarseGrid->getConverter().getLatticeTime(maxPhysT); ++iT) {
if (converge.hasConverged()) {
@@ -603,18 +676,18 @@ int main(int argc, char* argv[])
break;
}
- storeC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq);
+ fineCoupler.store();
coarseGrid->getSuperLattice().collideAndStream();
- interpolateStoredC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq);
fineGrid->getSuperLattice().collideAndStream();
- coupleC2F(*coarseGrid, *fineGrid, c2f_rho, c2f_u, c2f_fneq);
+ fineCoupler.interpolate();
+ fineCoupler.couple();
- storeC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq);
fineGrid->getSuperLattice().collideAndStream();
- coupleC2F(*coarseGrid, *fineGrid, c2f_rho, c2f_u, c2f_fneq);
+ fineCoupler.store();
+ fineCoupler.couple();
- coupleF2C(*coarseGrid, *fineGrid);
+ coarseCoupler.couple();
getResults(
"coarse_",