From 68f2384f79bff6553d1db21ac0b3173d57e4e2bf Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Thu, 10 Jan 2019 15:13:23 +0100 Subject: Add hacky MPI support for grid refinement Works but is nowhere near anything one could consider good. Obvious issues: * More than one cuboid per grid makes it harder to determine the next lattice cell to be coupled * i.e. currently lattice positions are determined ad hoc by resolving their physical position * Coupling is not actually parallelized * All coupling lines are traversed by all processes, way to much communication * Load balancing and cuboid decomposition doesn't care about refinement * ideally refined cuboids should be computationally near their coarse _parent_ cuboids The first two isses should be fixable with a reasonable amount of work. This sadly doesn't apply in any form to the last issue. --- apps/adrian/poiseuille2d/poiseuille2d.cpp | 13 +++--- src/refinement/coupler2D.hh | 77 +++++++++++++++++++++++-------- src/refinement/grid2D.hh | 7 ++- 3 files changed, 72 insertions(+), 25 deletions(-) diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp index a6d7bd1..df166e1 100644 --- a/apps/adrian/poiseuille2d/poiseuille2d.cpp +++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp @@ -35,10 +35,11 @@ typedef double T; #define DESCRIPTOR D2Q9Descriptor -const T lx = 4.0; // length of the channel -const T ly = 1.0; // height of the channel -const int N = 60; // resolution of the model +const T lx = 8.0; // length of the channel +const T ly = 2.0; // height of the channel +const int N = 64; // resolution of the model const T Re = 500.; // Reynolds number +const T baseTau = 0.57; // Relaxation time of coarsest grid 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 @@ -183,10 +184,10 @@ int main(int argc, char* argv[]) const Vector coarseExtend {lx, ly}; IndicatorCuboid2D coarseCuboid(coarseExtend, coarseOrigin); - auto coarseGrid = Grid2D::make(coarseCuboid, N, 0.57, Re); + auto coarseGrid = Grid2D::make(coarseCuboid, N, baseTau, Re); prepareGeometry(coarseGrid->getConverter(), coarseGrid->getSuperGeometry()); - const Vector wantedFineExtend {2.0, 0.75}; + const Vector wantedFineExtend {3.0, 1.5}; const Vector fineOrigin = coarseGrid->alignToGrid({0.8, (ly-wantedFineExtend[1])/2}); const Vector fineExtend = coarseGrid->alignToGrid(fineOrigin + wantedFineExtend) - fineOrigin; @@ -220,7 +221,7 @@ int main(int argc, char* argv[]) sOnLatticeBoundaryCondition2D fineSBoundaryCondition(fineGrid->getSuperLattice()); createLocalBoundaryCondition2D(fineSBoundaryCondition); - const Vector wantedFineExtend2 {0.4, 0.4}; + const Vector wantedFineExtend2 {0.6, 0.4}; const Vector fineOrigin2 = fineGrid->alignToGrid({1.05, (ly-wantedFineExtend2[1])/2}); const Vector fineExtend2 = fineGrid->alignToGrid(fineOrigin2 + wantedFineExtend2) - fineOrigin2; diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh index 8370707..ab082f4 100644 --- a/src/refinement/coupler2D.hh +++ b/src/refinement/coupler2D.hh @@ -35,10 +35,16 @@ template class DESCRIPTOR> Vector Coupler2D::getFineLatticeR(int y) const { if (_vertical) { - return _fineOrigin + Vector {0, 0, y}; + const auto physR = _physOrigin + Vector{0, y*_fine.getConverter().getPhysDeltaX()}; + Vector latticeR{}; + _fine.getCuboidGeometry().getLatticeR(physR, latticeR); + return latticeR; } else { - return _fineOrigin + Vector {0, y, 0}; + const auto physR = _physOrigin + Vector{y*_fine.getConverter().getPhysDeltaX(), 0}; + Vector latticeR{}; + _fine.getCuboidGeometry().getLatticeR(physR, latticeR); + return latticeR; } } @@ -46,10 +52,16 @@ template class DESCRIPTOR> Vector Coupler2D::getCoarseLatticeR(int y) const { if (_vertical) { - return _coarseOrigin + Vector {0, 0, y}; + const auto physR = _physOrigin + Vector{0, y*_coarse.getConverter().getPhysDeltaX()}; + Vector latticeR{}; + _coarse.getCuboidGeometry().getLatticeR(physR, latticeR); + return latticeR; } else { - return _coarseOrigin + Vector {0, y, 0}; + const auto physR = _physOrigin + Vector{y*_coarse.getConverter().getPhysDeltaX(), 0}; + Vector latticeR{}; + _coarse.getCuboidGeometry().getLatticeR(physR, latticeR); + return latticeR; } } @@ -97,12 +109,13 @@ void FineCoupler2D::store() 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); - coarseLattice.get(pos).computeFneq(fNeq); + Cell coarseCell; + coarseLattice.get(pos, coarseCell); + lbHelpers::computeRhoU(coarseCell, rho, u); + lbHelpers::computeFneq(coarseCell, fNeq, rho, u); _c2f_rho[y] = rho; _c2f_u[y] = Vector(u); @@ -134,17 +147,18 @@ void FineCoupler2D::interpolate() auto& coarseLattice = this->_coarse.getSuperLattice(); for (int y=0; y < this->_coarseSize; ++y) { - const auto coarseCell = coarseLattice.get(this->getCoarseLatticeR(y)); + Cell coarseCell; + coarseLattice.get(this->getCoarseLatticeR(y), coarseCell); T rho{}; T u[2] {}; - coarseCell.computeRhoU(rho, u); + 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]); T fNeq[DESCRIPTOR::q] {}; - coarseCell.computeFneq(fNeq); + 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]); } @@ -162,11 +176,16 @@ void FineCoupler2D::couple() const auto finePos = this->getFineLatticeR(2*y); T fEq[DESCRIPTOR::q] {}; - coarseLattice.get(coarsePos).computeFeq(fEq); + Cell coarseCell; + coarseLattice.get(coarsePos, coarseCell); + lbHelpers::computeFeq(coarseCell, fEq); + Cell cell; + fineLattice.get(finePos, cell); for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - fineLattice.get(finePos)[iPop] = fEq[iPop] + this->_coarse.getScalingFactor() * _c2f_fneq[y][iPop]; + cell[iPop] = fEq[iPop] + this->_coarse.getScalingFactor() * _c2f_fneq[y][iPop]; } + fineLattice.set(finePos, cell); } for (int y=1; y < this->_coarseSize-2; ++y) { @@ -192,6 +211,8 @@ void FineCoupler2D::couple() const T uSqr = u[0]*u[0] + u[1]*u[1]; 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( @@ -201,8 +222,10 @@ void FineCoupler2D::couple() _c2f_fneq[y+2][iPop] ); - fineLattice.get(finePos)[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq; + fineCell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq; } + + fineLattice.set(finePos, fineCell); } { @@ -225,6 +248,8 @@ void FineCoupler2D::couple() const T uSqr = u[0]*u[0] + u[1]*u[1]; const auto finePos = this->getFineLatticeR(1); + Cell fineCell; + fineLattice.get(finePos, fineCell); for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { const T fneq = order3interpolation( @@ -232,8 +257,10 @@ void FineCoupler2D::couple() _c2f_fneq[1][iPop], _c2f_fneq[2][iPop] ); - fineLattice.get(finePos)[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq; + fineCell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq; } + + fineLattice.set(finePos, fineCell); } { @@ -256,6 +283,8 @@ void FineCoupler2D::couple() const T uSqr = u[0]*u[0] + u[1]*u[1]; 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( @@ -263,8 +292,10 @@ void FineCoupler2D::couple() _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; + fineCell[iPop] = lbHelpers::equilibrium(iPop, rho, u, uSqr) + this->_coarse.getScalingFactor() * fneq; } + + fineLattice.set(finePos, fineCell); } } @@ -276,7 +307,10 @@ void computeRestrictedFneq(const SuperLattice2D& lattice, { for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { T fNeq[DESCRIPTOR::q] {}; - lattice.get(pos[0], pos[1] + DESCRIPTOR::c[iPop][0], pos[2] + DESCRIPTOR::c[iPop][1]).computeFneq(fNeq); + Cell cell; + const Vector neighborPos {pos[0], pos[1] + DESCRIPTOR::c[iPop][0], pos[2] + DESCRIPTOR::c[iPop][1]}; + lattice.get(neighborPos, cell); + lbHelpers::computeFneq(cell, fNeq); for (int jPop=0; jPop < DESCRIPTOR::q; ++jPop) { restrictedFneq[jPop] += fNeq[jPop]; @@ -311,14 +345,21 @@ void CoarseCoupler2D::couple() const auto coarsePos = this->getCoarseLatticeR(y); T fEq[DESCRIPTOR::q] {}; - fineLattice.get(finePos).computeFeq(fEq); + Cell fineCell; + fineLattice.get(finePos, fineCell); + lbHelpers::computeFeq(fineCell, fEq); T fNeq[DESCRIPTOR::q] {}; computeRestrictedFneq(fineLattice, finePos, fNeq); + Cell coarseCell; + coarseLattice.get(coarsePos, coarseCell); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - coarseLattice.get(coarsePos)[iPop] = fEq[iPop] + this->_coarse.getInvScalingFactor() * fNeq[iPop]; + coarseCell[iPop] = fEq[iPop] + this->_coarse.getInvScalingFactor() * fNeq[iPop]; } + + coarseLattice.set(coarsePos, coarseCell); } } diff --git a/src/refinement/grid2D.hh b/src/refinement/grid2D.hh index 71b943e..65dd84d 100644 --- a/src/refinement/grid2D.hh +++ b/src/refinement/grid2D.hh @@ -55,7 +55,12 @@ Grid2D::Grid2D(IndicatorF2D& domainF, int resolution, T tau, in _cuboids(new CuboidGeometry2D( domainF, _converter->getConversionFactorLength(), - 1)), +#ifdef PARALLEL_MODE_MPI + singleton::mpi().getSize() +#else + 1 +#endif + )), _balancer(new HeuristicLoadBalancer( *_cuboids)), _geometry(new SuperGeometry2D( -- cgit v1.2.3