From 78d8d5206b0986e81690d2cee6b7202abef5a1f2 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Tue, 8 Jan 2019 10:43:28 +0100 Subject: Extract refinement scaffolding into separate units --- apps/adrian/poiseuille2d/poiseuille2d.cpp | 496 +----------------------------- 1 file changed, 1 insertion(+), 495 deletions(-) (limited to 'apps/adrian/poiseuille2d') diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp index 9ef1ce8..f30e9b4 100644 --- a/apps/adrian/poiseuille2d/poiseuille2d.cpp +++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp @@ -197,500 +197,6 @@ void getResults(const std::string& prefix, } } -template class DESCRIPTOR> class FineCoupler; -template class DESCRIPTOR> class CoarseCoupler; - -template class DESCRIPTOR> -class Grid { -private: - std::unique_ptr> _converter; - std::unique_ptr> _cuboids; - std::unique_ptr> _balancer; - std::unique_ptr> _geometry; - std::unique_ptr> _lattice; - - std::vector>> _fineGrids; - - std::vector>> _fineCouplers; - std::vector>> _coarseCouplers; - -public: - static std::unique_ptr> make(IndicatorF2D& domainF, int resolution, T tau, int re) - { - return std::unique_ptr>( - new Grid(domainF, resolution, tau, re) - ); - } - - Grid(IndicatorF2D& domainF, int resolution, T tau, int re): - _converter(new UnitConverterFromResolutionAndRelaxationTime( - resolution, // resolution: number of voxels per charPhysL - tau, // latticeRelaxationTime: relaxation time, has to be greater than 0.5! - T{1}, // charPhysLength: reference length of simulation geometry - T{1}, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__ - T{1./re}, // physViscosity: physical kinematic viscosity in __m^2 / s__ - T{1}, // physDensity: physical density in __kg / m^3__ - T{1})), - _cuboids(new CuboidGeometry2D( - domainF, - _converter->getConversionFactorLength(), - 1)), - _balancer(new HeuristicLoadBalancer( - *_cuboids)), - _geometry(new SuperGeometry2D( - *_cuboids, - *_balancer, - 2)), - _lattice(new SuperLattice2D( - *_geometry)) - { - _converter->print(); - } - - UnitConverter& getConverter() - { - return *_converter; - } - - CuboidGeometry2D& getCuboidGeometry() - { - return *_cuboids; - } - - LoadBalancer& getLoadBalancer() - { - return *_balancer; - } - - SuperGeometry2D& getSuperGeometry() - { - return *_geometry; - } - - SuperLattice2D& getSuperLattice() - { - return *_lattice; - } - - T getScalingFactor() - { - return 4. - _converter->getLatticeRelaxationFrequency(); - } - - T getInvScalingFactor() - { - return 1./getScalingFactor(); - } - - void collideAndStream() - { - for ( auto& fineCoupler : _fineCouplers ) { - fineCoupler->store(); - } - - this->getSuperLattice().collideAndStream(); - - for ( auto& fineGrid : _fineGrids ) { - fineGrid->getSuperLattice().collideAndStream(); - } - - for ( auto& fineCoupler : _fineCouplers ) { - fineCoupler->interpolate(); - fineCoupler->couple(); - } - - for ( auto& fineGrid : _fineGrids ) { - fineGrid->getSuperLattice().collideAndStream(); - } - - for ( auto& fineCoupler : _fineCouplers ) { - fineCoupler->store(); - fineCoupler->couple(); - } - - for ( auto& coarseCoupler : _coarseCouplers ) { - coarseCoupler->couple(); - } - } - - FineCoupler& addFineCoupling( - Grid& fineGrid, Vector origin, Vector extend) - { - _fineCouplers.emplace_back( - new FineCoupler( - *this, fineGrid, origin, extend)); - return *_fineCouplers.back(); - } - - CoarseCoupler& addCoarseCoupling( - Grid& fineGrid, Vector origin, Vector extend) - { - _coarseCouplers.emplace_back( - new CoarseCoupler( - *this, fineGrid, origin, extend)); - return *_coarseCouplers.back(); - } - - Grid& refine(IndicatorF2D& domainF) - { - _fineGrids.emplace_back( - new Grid( - domainF, - 2*getConverter().getResolution(), - 2.0*getConverter().getLatticeRelaxationTime() - 0.5, - getConverter().getReynoldsNumber() - )); - return *_fineGrids.back(); - } - - Grid& refine(Vector origin, Vector extend) - { - IndicatorCuboid2D fineCuboid(extend, origin); - auto& fineGrid = refine(fineCuboid); - - const Vector extendX = {extend[0],0}; - const Vector extendY = {0,extend[1]}; - - addFineCoupling(fineGrid, origin, extendY); - addFineCoupling(fineGrid, origin + extendX, extendY); - addFineCoupling(fineGrid, origin + extendY, extendX); - addFineCoupling(fineGrid, origin, extendX); - - const T coarseDeltaX = getConverter().getPhysDeltaX(); - const Vector innerOrigin = origin + coarseDeltaX; - const Vector innerExtendX = extendX - Vector {2*coarseDeltaX,0}; - const Vector innerExtendY = extendY - Vector {0,2*coarseDeltaX}; - - addCoarseCoupling(fineGrid, innerOrigin, innerExtendY); - addCoarseCoupling(fineGrid, innerOrigin + innerExtendX, innerExtendY); - addCoarseCoupling(fineGrid, innerOrigin + innerExtendY, innerExtendX); - addCoarseCoupling(fineGrid, innerOrigin, innerExtendX); - - return fineGrid; - } -}; - -template class DESCRIPTOR> -void computeFeq(const Cell& cell, T fEq[DESCRIPTOR::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::q; ++iPop) { - fEq[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr); - } -} - -template class DESCRIPTOR> -void computeFneq(const Cell& cell, T fNeq[DESCRIPTOR::q]) -{ - T rho{}; - T u[2] {}; - cell.computeRhoU(rho, u); - lbHelpers::computeFneq(cell, fNeq, rho, u); -} - -T order4interpolation(T fm3, T fm1, T f1, T f3) -{ - return 9./16. * (fm1 + f1) - 1./16. * (fm3 + f3); -} - -T order3interpolation(T fm1, T f1, T f3) -{ - return 3./8. * fm1 + 3./4. * f1 - 1./8. * f3; -} - -T order2interpolation(T f0, T f1) -{ - return 0.5 * (f0 + f1); -} - -template class DESCRIPTOR> -class Coupler { -protected: - Grid& _coarse; - Grid& _fine; - - const int _coarseSize; - const int _fineSize; - const bool _vertical; - - Vector _physOrigin; - Vector _coarseOrigin; - Vector _fineOrigin; - - Vector getFineLatticeR(int y) const - { - if (_vertical) { - return _fineOrigin + Vector {0, 0, y}; - } - else { - return _fineOrigin + Vector {0, y, 0}; - } - } - - Vector getCoarseLatticeR(int y) const - { - if (_vertical) { - return _coarseOrigin + Vector {0, 0, y}; - } - else { - return _coarseOrigin + Vector {0, y, 0}; - } - } - -public: - Coupler(Grid& coarse, Grid& fine, Vector origin, Vector extend): - _coarse(coarse), - _fine(fine), - _coarseSize(floor(extend.norm() / coarse.getConverter().getPhysDeltaX() + 0.5)+1), - _fineSize(2*_coarseSize-1), - _vertical(util::nearZero(extend[0])) - { - OLB_ASSERT(util::nearZero(extend[0]) || util::nearZero(extend[1]), "Coupling is only implemented alongside unit vectors"); - - const auto& coarseGeometry = _coarse.getCuboidGeometry(); - const auto& fineGeometry = _fine.getCuboidGeometry(); - - Vector tmpLatticeOrigin{}; - coarseGeometry.getLatticeR(origin, tmpLatticeOrigin); - _physOrigin = coarseGeometry.getPhysR(tmpLatticeOrigin.toStdVector()); - - coarseGeometry.getLatticeR(_physOrigin, _coarseOrigin); - fineGeometry.getLatticeR(_physOrigin, _fineOrigin); - } -}; - -template class DESCRIPTOR> -class FineCoupler : public Coupler { -private: - std::vector _c2f_rho; - std::vector> _c2f_u; - std::vector::q>> _c2f_fneq; - -public: - FineCoupler(Grid& coarse, Grid& fine, Vector origin, Vector extend): - Coupler(coarse, fine, origin, extend), - _c2f_rho(this->_coarseSize), - _c2f_u(this->_coarseSize, Vector(T{})), - _c2f_fneq(this->_coarseSize, Vector::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: " << this->_fineSize << std::endl; - } - - void store() - { - auto& coarseLattice = this->_coarse.getSuperLattice(); - - for (int y=0; y < this->_coarseSize; ++y) { - const auto pos = this->getCoarseLatticeR(y); - - T rho{}; - T u[2] {}; - T fNeq[DESCRIPTOR::q] {}; - coarseLattice.get(pos).computeRhoU(rho, u); - computeFneq(coarseLattice.get(pos), fNeq); - - _c2f_rho[y] = rho; - _c2f_u[y] = Vector(u); - _c2f_fneq[y] = Vector::q>(fNeq); - } - } - - void interpolate() - { - 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::q] {}; - computeFneq(coarseLattice.get(coarsePos), fNeq); - for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - _c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], _c2f_fneq[y][iPop]); - } - } - } - - void couple() - { - 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::q] {}; - computeFeq(coarseLattice.get(coarsePos), fEq); - - for (int iPop=0; iPop < DESCRIPTOR::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::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::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::q; ++iPop) { - const T fneq = order3interpolation( - _c2f_fneq[0][iPop], - _c2f_fneq[1][iPop], - _c2f_fneq[2][iPop] - ); - fineLattice.get(finePos)[iPop] = lbHelpers::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(this->_fineSize-2); - - 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] - ); - fineLattice.get(finePos)[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq; - } - } - } -}; - -template class DESCRIPTOR> -void computeRestrictedFneq(const SuperLattice2D& lattice, Vector pos, T restrictedFneq[DESCRIPTOR::q]) -{ - for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - T fNeq[DESCRIPTOR::q] {}; - computeFneq(lattice.get(pos[0], pos[1] + DESCRIPTOR::c[iPop][0], pos[2] + DESCRIPTOR::c[iPop][1]), fNeq); - - for (int jPop=0; jPop < DESCRIPTOR::q; ++jPop) { - restrictedFneq[jPop] += fNeq[jPop]; - } - } - - for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - restrictedFneq[iPop] /= DESCRIPTOR::q; - } -} - -template class DESCRIPTOR> -class CoarseCoupler : public Coupler { -public: - CoarseCoupler(Grid& coarse, Grid& fine, Vector origin, Vector extend): - Coupler(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(); - - for (int y=0; y < this->_coarseSize; ++y) { - const auto finePos = this->getFineLatticeR(2*y); - const auto coarsePos = this->getCoarseLatticeR(y); - - T fEq[DESCRIPTOR::q] {}; - computeFeq(fineLattice.get(finePos), fEq); - - T fNeq[DESCRIPTOR::q] {}; - computeRestrictedFneq(fineLattice, finePos, fNeq); - - for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - coarseLattice.get(coarsePos)[iPop] = fEq[iPop] + this->_coarse.getInvScalingFactor() * fNeq[iPop]; - } - } - } -}; - int main(int argc, char* argv[]) { olbInit(&argc, &argv); @@ -704,7 +210,7 @@ int main(int argc, char* argv[]) Vector fineOrigin { 0.25 * lx, 0.2 }; Vector fineExtend { 0.6 * lx, 0.6 }; - auto coarseGrid = Grid::make(coarseCuboid, N, 0.8, Re); + auto coarseGrid = Grid2D::make(coarseCuboid, N, 0.8, Re); auto fineGrid = &coarseGrid->refine(fineOrigin, fineExtend); prepareGeometry(coarseGrid->getConverter(), coarseGrid->getSuperGeometry()); -- cgit v1.2.3