From 9887f898074efb6a2682eedb174b77d9b9b63e0c Mon Sep 17 00:00:00 2001 From: Rupert Nash Date: Wed, 31 Jan 2024 15:03:16 +0000 Subject: [PATCH] Switch domain and LBM to use indexes not counts for sites of different collision types --- Code/geometry/Domain.cc | 214 +++------ Code/geometry/Domain.h | 91 ++-- Code/lb/CMakeLists.txt | 1 + Code/lb/EntropyTester.h | 297 +++++------- Code/lb/LbmParameters.h | 3 - Code/lb/lb.cc | 73 +++ Code/lb/lb.h | 157 ++++--- Code/lb/lb.hpp | 385 +++++---------- Code/lb/streamers/BulkStreamer.h | 6 +- Code/lb/streamers/JunkYang.h | 8 +- Code/lb/streamers/StreamerTypeFactory.h | 8 +- Code/lb/streamers/VirtualSiteIolet.h | 8 +- Code/multiscale/MultiscaleSimulationMaster.h | 61 +-- .../tests/helpers/FourCubeBasedTestFixture.cc | 1 - Code/tests/lb/GuoForcingTests.cc | 8 +- Code/tests/lb/StreamerTests.cc | 440 ++++++++---------- .../tests/lb/VirtualSiteIoletStreamerTests.cc | 70 +-- 17 files changed, 791 insertions(+), 1040 deletions(-) create mode 100644 Code/lb/lb.cc diff --git a/Code/geometry/Domain.cc b/Code/geometry/Domain.cc index bc7f36c19..ce9351843 100644 --- a/Code/geometry/Domain.cc +++ b/Code/geometry/Domain.cc @@ -18,14 +18,14 @@ namespace hemelb::geometry Domain::Domain(const lb::LatticeInfo& latticeInfo, const net::IOCommunicator& comms_) : latticeInfo(latticeInfo), - shared_counts(comms_, 0), + shared_site_type_indices(comms_, 0), neighbouringData(new neighbouring::NeighbouringDomain(latticeInfo)), comms(comms_) { } Domain::Domain(const lb::LatticeInfo& latticeInfo, GmyReadResult& readResult, const net::IOCommunicator& comms_) : - latticeInfo(latticeInfo), shared_counts(comms_, 0), + latticeInfo(latticeInfo), shared_site_type_indices(comms_, 0), neighbouringData(new neighbouring::NeighbouringDomain(latticeInfo)), rank_for_site_store(std::move(readResult.block_store)), comms(comms_) @@ -69,21 +69,14 @@ namespace hemelb::geometry void Domain::ProcessReadSites(const GmyReadResult & readResult) { log::Logger::Log("Processing sites assigned to each MPI process"); + const auto NO_NORMAL = util::Vector3D(NO_VALUE); const auto max_site_index = sites - util::Vector3D::Ones(); blocks.resize(rank_for_site_store->GetBlockCount()); totalSharedFs = 0; - std::vector domainEdgeSiteData[COLLISION_TYPES]; - std::vector midDomainSiteData[COLLISION_TYPES]; - std::vector domainEdgeBlockNumber[COLLISION_TYPES]; - std::vector midDomainBlockNumber[COLLISION_TYPES]; - std::vector domainEdgeSiteNumber[COLLISION_TYPES]; - std::vector midDomainSiteNumber[COLLISION_TYPES]; - std::vector > domainEdgeWallNormals[COLLISION_TYPES]; - std::vector > midDomainWallNormals[COLLISION_TYPES]; - std::vector domainEdgeWallDistance[COLLISION_TYPES]; - std::vector midDomainWallDistance[COLLISION_TYPES]; + ReadDataBundle domainEdge; + ReadDataBundle midDomain; proc_t localRank = comms.Rank(); @@ -153,7 +146,8 @@ namespace hemelb::geometry if (blocks[blockOctIdx].IsEmpty()) blocks[blockOctIdx] = Block(GetSitesPerBlockVolumeUnit()); - auto assignedRank = blockReadIn.Sites[localSiteId].targetProcessor; + auto const& siteReadIn = blockReadIn.Sites[localSiteId]; + auto assignedRank = siteReadIn.targetProcessor; blocks[blockOctIdx].SetProcessorRankForSite(localSiteId, assignedRank); // If the site is not on this processor, continue. @@ -164,7 +158,7 @@ namespace hemelb::geometry util::Vector3D siteGlobalCoords = lowest_site_in_block + siteTraverser.GetCurrentLocation(); // Iterate over all non-zero direction vectors. for (unsigned int l = 1; l < latticeInfo.GetNumVectors(); l++) { - // Find the neighbour site co-ords in this direction. + // Find the neighbour site coords in this direction. util::Vector3D neighbourGlobalCoords = siteGlobalCoords + latticeInfo.GetVector(l).as(); Vec16 neighbourBlock; @@ -183,7 +177,7 @@ namespace hemelb::geometry // Set the collision type data. map_block site data is renumbered according to // fluid site numbers within a particular collision type. - SiteData siteData(blockReadIn.Sites[localSiteId]); + SiteData siteData(siteReadIn); int l = -1; switch (siteData.GetCollisionType()) { case FLUID: @@ -206,40 +200,26 @@ namespace hemelb::geometry break; } - const util::Vector3D& normal = blockReadIn.Sites[localSiteId].wallNormalAvailable ? - blockReadIn.Sites[localSiteId].wallNormal : - util::Vector3D(NO_VALUE); + auto& normal = siteReadIn.wallNormalAvailable ? siteReadIn.wallNormal : NO_NORMAL; - if (isMidDomainSite) { - midDomainBlockNumber[l].push_back(blockOctIdx); - midDomainSiteNumber[l].push_back(localSiteId); - midDomainSiteData[l].push_back(siteData); - midDomainWallNormals[l].push_back(normal); + auto push_site_onto_bundle = [&](ReadDataBundle& b) { + b.blockNumbers[l].push_back(blockOctIdx); + b.siteNumbers[l].push_back(localSiteId); + b.siteData[l].push_back(siteData); + b.wallNormals[l].push_back(normal); for (Direction direction = 1; direction < latticeInfo.GetNumVectors(); direction++) { - midDomainWallDistance[l].push_back(blockReadIn.Sites[localSiteId].links[direction - 1].distanceToIntersection); + b.wallDistance[l].push_back(siteReadIn.links[direction - 1].distanceToIntersection); } + }; + if (isMidDomainSite) { + push_site_onto_bundle(midDomain); } else { - domainEdgeBlockNumber[l].push_back(blockOctIdx); - domainEdgeSiteNumber[l].push_back(localSiteId); - domainEdgeSiteData[l].push_back(siteData); - domainEdgeWallNormals[l].push_back(normal); - for (Direction direction = 1; direction < latticeInfo.GetNumVectors(); direction++) { - domainEdgeWallDistance[l].push_back(blockReadIn.Sites[localSiteId].links[direction - 1].distanceToIntersection); - } + push_site_onto_bundle(domainEdge); } } } - PopulateWithReadData(midDomainBlockNumber, - midDomainSiteNumber, - midDomainSiteData, - midDomainWallNormals, - midDomainWallDistance, - domainEdgeBlockNumber, - domainEdgeSiteNumber, - domainEdgeSiteData, - domainEdgeWallNormals, - domainEdgeWallDistance); + PopulateWithReadData(midDomain, domainEdge); // We now have the distributed store setup, so can find which // process any site lives one. Set up the neighbouring @@ -285,76 +265,51 @@ namespace hemelb::geometry } } - void Domain::PopulateWithReadData( - const std::vector midDomainBlockNumbers[COLLISION_TYPES], - const std::vector midDomainSiteNumbers[COLLISION_TYPES], - const std::vector midDomainSiteData[COLLISION_TYPES], - const std::vector > midDomainWallNormals[COLLISION_TYPES], - const std::vector midDomainWallDistance[COLLISION_TYPES], - const std::vector domainEdgeBlockNumbers[COLLISION_TYPES], - const std::vector domainEdgeSiteNumbers[COLLISION_TYPES], - const std::vector domainEdgeSiteData[COLLISION_TYPES], - const std::vector > domainEdgeWallNormals[COLLISION_TYPES], - const std::vector domainEdgeWallDistance[COLLISION_TYPES]) - { - auto write_my_sites = rank_for_site_store->begin_writes(); - - log::Logger::Log("Assigning local indices to sites and associated data"); - // Populate the collision count arrays. - for (unsigned collisionType = 0; collisionType < COLLISION_TYPES; collisionType++) - { - MidDomainCollisionCount(collisionType) = midDomainBlockNumbers[collisionType].size(); - DomainEdgeCollisionCount(collisionType) = domainEdgeBlockNumbers[collisionType].size(); - } - // Data about local sites. - SiteRankIndex rank_index = {comms.Rank(), 0}; + void Domain::PopulateWithReadData(const ReadDataBundle &midDomain, const ReadDataBundle &domainEdge) + { + auto write_my_sites = rank_for_site_store->begin_writes(); + + log::Logger::Log("Assigning local indices to sites and associated data"); + + // Data about local sites. + SiteRankIndex rank_index = {comms.Rank(), 0}; + // Track index of changes in type and store them in the counts array + auto ranges = shared_site_type_indices.Span(); + auto current_count_iter = begin(ranges); + // Merge local sites into one big array. This lambda does either all the + // midDomain or all the domainEdge sites, putting them in order of their + // collision type (bulk/wall/inlet/etc). + auto op = [this, &rank_index, &write_my_sites, ¤t_count_iter](ReadDataBundle const& b) { + auto const NV = latticeInfo.GetNumVectors(); auto& localFluidSites = rank_index[1]; - // Data about contiguous local sites. First midDomain stuff, then domainEdge. - for (unsigned collisionType = 0; collisionType < COLLISION_TYPES; collisionType++) - { - for (unsigned indexInType = 0; indexInType < GetMidDomainCollisionCount(collisionType); - indexInType++) - { - siteData.push_back(midDomainSiteData[collisionType][indexInType]); - wallNormalAtSite.emplace_back(midDomainWallNormals[collisionType][indexInType]); - for (Direction direction = 1; direction < latticeInfo.GetNumVectors(); direction++) - { - distanceToWall.push_back(midDomainWallDistance[collisionType][indexInType - * (latticeInfo.GetNumVectors() - 1) + direction - 1]); + for (unsigned collisionType = 0; collisionType < COLLISION_TYPES; collisionType++) { + *current_count_iter = localFluidSites; + + for (unsigned i = 0; i < b.blockNumbers[collisionType].size(); ++i) { + siteData.push_back(b.siteData[collisionType][i]); + wallNormalAtSite.emplace_back(b.wallNormals[collisionType][i]); + for (Direction dir = 1; dir < NV; dir++) { + distanceToWall.push_back( + b.wallDistance[collisionType][i * (NV - 1) + dir - 1] + ); } - site_t blockId = midDomainBlockNumbers[collisionType][indexInType]; - site_t siteId = midDomainSiteNumbers[collisionType][indexInType]; + site_t blockId = b.blockNumbers[collisionType][i]; + site_t siteId = b.siteNumbers[collisionType][i]; blocks[blockId].SetLocalContiguousIndexForSite(siteId, localFluidSites); globalSiteCoords.push_back(GetGlobalCoords(blockId, GetSiteCoordsFromSiteId(siteId))); write_my_sites(blockId)(siteId) = rank_index; localFluidSites++; } + // Note for all iterations except the last we will overwrite this, but nice and clear. + ++current_count_iter; + *current_count_iter = localFluidSites; } - - for (unsigned collisionType = 0; collisionType < COLLISION_TYPES; collisionType++) - { - for (unsigned indexInType = 0; indexInType < GetDomainEdgeCollisionCount(collisionType); - indexInType++) - { - siteData.push_back(domainEdgeSiteData[collisionType][indexInType]); - wallNormalAtSite.emplace_back(domainEdgeWallNormals[collisionType][indexInType]); - for (Direction direction = 1; direction < latticeInfo.GetNumVectors(); direction++) - { - distanceToWall.push_back(domainEdgeWallDistance[collisionType][indexInType - * (latticeInfo.GetNumVectors() - 1) + direction - 1]); - } - site_t blockId = domainEdgeBlockNumbers[collisionType][indexInType]; - site_t siteId = domainEdgeSiteNumbers[collisionType][indexInType]; - blocks[blockId].SetLocalContiguousIndexForSite(siteId, localFluidSites); - globalSiteCoords.push_back(GetGlobalCoords(blockId, GetSiteCoordsFromSiteId(siteId))); - write_my_sites(blockId)(siteId) = rank_index; - localFluidSites++; - } - - } - LocalFluidSiteCount() = localFluidSites; - } + }; + // Now apply it. Order of these is important! + op(midDomain); + op(domainEdge); + } void Domain::CollectFluidSiteDistribution() { @@ -362,7 +317,7 @@ namespace hemelb::geometry fluidSitesOnEachProcessor = comms.AllGather(GetLocalFluidSiteCount()); totalFluidSites = std::reduce( fluidSitesOnEachProcessor.begin(), fluidSitesOnEachProcessor.end(), - 0, std::plus<>{}); + site_t(0), std::plus<>{}); } void Domain::CollectGlobalSiteExtrema() @@ -691,64 +646,23 @@ namespace hemelb::geometry siteCoords = (location % blockSize).as(); } - template - site_t total(Iter&& begin) { - return std::reduce(begin, begin + COLLISION_TYPES, site_t(0), std::plus<>()); - } - bool Domain::IsSiteDomainEdge(int rank, site_t local_idx) const { if (!remote_counts_cache.contains(rank)) { auto &tmp = remote_counts_cache[rank]; - shared_counts.Get(std::span{tmp}, rank); + shared_site_type_indices.Get(std::span{tmp}, rank); } auto const& counts = remote_counts_cache[rank]; - // Mid-domain counts are first - auto n_mid_domain = total(counts.begin()); + // Mid-domain indices are first + auto n_mid_domain = counts[COLLISION_TYPES]; return local_idx >= n_mid_domain; } - site_t Domain::GetMidDomainSiteCount() const + Vec16 Domain::GetBlockIJK(site_t block) const { - return total(&GetMidDomainCollisionCount(0)); - } - site_t Domain::GetDomainEdgeSiteCount() const { - return total(&GetDomainEdgeCollisionCount(0)); + auto& t = rank_for_site_store->GetTree(); + return t.GetLeafCoords(block); } - Vec16 Domain::GetBlockIJK(site_t block) const - { - auto& t = rank_for_site_store->GetTree(); - return t.GetLeafCoords(block); - } - -/* void FieldData::SendAndReceive(hemelb::net::Net* net) - { - for (std::vector::const_iterator it = neighbouringProcs.begin(); - it != neighbouringProcs.end(); ++it) - { - // Request the receive into the appropriate bit of FOld. - net->RequestReceive(GetFOld( (*it).FirstSharedDistribution), - (int) ( ( (*it).SharedDistributionCount)), - (*it).Rank); - // Request the send from the right bit of FNew. - net->RequestSend(GetFNew( (*it).FirstSharedDistribution), - (int) ( ( (*it).SharedDistributionCount)), - (*it).Rank); - - } - }*/ - - /*void FieldData::CopyReceived() - { - // Copy the distribution functions received from the neighbouring - // processors into the destination buffer "f_new". - for (site_t i = 0; i < totalSharedFs; i++) - { - *GetFNew(streamingIndicesForReceivedDistributions[i]) = - *GetFOld(neighbouringProcs[0].FirstSharedDistribution + i); - } - }*/ - void Domain::Report(reporting::Dict& dictionary) { dictionary.SetIntValue("SITES", GetTotalFluidSites()); diff --git a/Code/geometry/Domain.h b/Code/geometry/Domain.h index 60780cace..70137fb49 100644 --- a/Code/geometry/Domain.h +++ b/Code/geometry/Domain.h @@ -58,6 +58,7 @@ namespace hemelb::geometry friend class tests::helpers::LatticeDataAccess; public: + template friend class lb::LBM; //! Let the LBM have access to internals so it can initialise the distribution arrays. template @@ -179,7 +180,7 @@ namespace hemelb::geometry */ inline site_t const& GetLocalFluidSiteCount() const { - return shared_counts.Span()[2 * COLLISION_TYPES]; + return shared_site_type_indices.Span()[2 * COLLISION_TYPES]; } site_t GetContiguousSiteId(util::Vector3D location) const; @@ -254,8 +255,30 @@ namespace hemelb::geometry // Will do RMA iff required data not cached. bool IsSiteDomainEdge(int rank, site_t local_idx) const; - site_t GetMidDomainSiteCount() const; - site_t GetDomainEdgeSiteCount() const; + inline site_t GetMidDomainSiteCount() const { + return shared_site_type_indices.Span()[COLLISION_TYPES]; + } + inline site_t GetDomainEdgeSiteCount() const { + return shared_site_type_indices.Span()[2 * COLLISION_TYPES] - GetMidDomainSiteCount(); + } + + inline std::pair GetMidDomainSiteRange(unsigned ct) const { + auto&& span = shared_site_type_indices.Span(); + return {span[ct], span[ct+1]}; + } + inline std::pair GetDomainEdgeSiteRange(unsigned ct) const { + auto&& span = shared_site_type_indices.Span(); + return {span[COLLISION_TYPES + ct], span[COLLISION_TYPES + ct + 1]}; + } + + inline std::pair GetAllMidDomainSiteRange() const { + auto&& span = shared_site_type_indices.Span(); + return {span[0], span[COLLISION_TYPES]}; + } + inline std::pair GetAllDomainEdgeSiteRange() const { + auto&& span = shared_site_type_indices.Span(); + return {span[COLLISION_TYPES], span[2 * COLLISION_TYPES]}; + } /** * Number of sites with all fluid neighbours residing on this rank, for the given @@ -263,9 +286,10 @@ namespace hemelb::geometry * @param collisionType * @return */ - inline site_t const& GetMidDomainCollisionCount(unsigned int collisionType) const + inline site_t GetMidDomainSiteCount(unsigned int collisionType) const { - return shared_counts.Span()[collisionType]; + auto [begin, end] = GetMidDomainSiteRange(collisionType); + return end - begin; } /** @@ -274,26 +298,17 @@ namespace hemelb::geometry * @param collisionType * @return */ - inline site_t const& GetDomainEdgeCollisionCount(unsigned int collisionType) const + inline site_t GetDomainEdgeSiteCount(unsigned int collisionType) const { - return shared_counts.Span()[COLLISION_TYPES + collisionType]; - } + auto [begin, end] = GetDomainEdgeSiteRange(collisionType); + return end - begin; - private: - // The setters should be private - inline site_t& MidDomainCollisionCount(unsigned int collisionType) - { - return shared_counts.Span()[collisionType]; - } - inline site_t& DomainEdgeCollisionCount(unsigned int collisionType) - { - return shared_counts.Span()[COLLISION_TYPES + collisionType]; } + inline site_t& LocalFluidSiteCount() { - return shared_counts.Span()[2 * COLLISION_TYPES]; + return shared_site_type_indices.Span()[2 * COLLISION_TYPES]; } - public: /** * Get the number of fluid sites on the given rank. @@ -342,7 +357,7 @@ namespace hemelb::geometry inline Site GetSite(site_t i) const { return Site{i, *this}; } - protected: + protected: /** * The protected default constructor does nothing. It exists to allow derivation from this * class for the purpose of testing. @@ -354,18 +369,21 @@ namespace hemelb::geometry void ProcessReadSites(const GmyReadResult& readResult); - // TODO: should these parameters be copies? + template + using ArrayOfVectors = std::array, COLLISION_TYPES>; + + struct ReadDataBundle { + ArrayOfVectors blockNumbers; + ArrayOfVectors siteNumbers; + ArrayOfVectors siteData; + ArrayOfVectors > wallNormals; + ArrayOfVectors wallDistance; + }; + void PopulateWithReadData( - const std::vector midDomainBlockNumbers[COLLISION_TYPES], - const std::vector midDomainSiteNumbers[COLLISION_TYPES], - const std::vector midDomainSiteData[COLLISION_TYPES], - const std::vector > midDomainWallNormals[COLLISION_TYPES], - const std::vector midDomainWallDistance[COLLISION_TYPES], - const std::vector domainEdgeBlockNumbers[COLLISION_TYPES], - const std::vector domainEdgeSiteNumbers[COLLISION_TYPES], - const std::vector domainEdgeSiteData[COLLISION_TYPES], - const std::vector > domainEdgeWallNormals[COLLISION_TYPES], - const std::vector domainEdgeWallDistance[COLLISION_TYPES]); + ReadDataBundle const& midDomain, + ReadDataBundle const& domainEdge + ); void CollectFluidSiteDistribution(); void CollectGlobalSiteExtrema(); @@ -503,11 +521,14 @@ namespace hemelb::geometry std::vector neighbouringProcs; //! Info about processors with neighbouring fluid sites. // In this window, we store the counts of stuff that need to be PGAS accessible. - // First COLLISION_TYPES midDomainProcCollisions; //! Number of fluid sites with all fluid neighbours on this rank, for each collision type. - // Second COLLISION_TYPES> domainEdgeProcCollisions; //! Number of fluid sites with at least one fluid neighbour on another rank, for each collision type. - // Last localFluidSites; //! The number of local fluid sites. + // A mid-domain site has all its neighbours also on this rank + // A domain edge site has at least one neighbour on another rank + // First, the start index of the mid-domain sites, grouped by their collision type. + // Next, the start index of the domain-edge sites, grouped by the collision type. + // Finally, the total number of local fluid sites. + // Thus index i + 1 holds the past-the-end index for range i static constexpr MPI_Aint N_SHARED = 2 * COLLISION_TYPES + 1; - net::WinData shared_counts; + net::WinData shared_site_type_indices; // Here we cache this data locally, since we expect to need // only a small number of other process's data, but many times using SharedCountArray = std::array; diff --git a/Code/lb/CMakeLists.txt b/Code/lb/CMakeLists.txt index 1dcc66817..93d56e6ea 100644 --- a/Code/lb/CMakeLists.txt +++ b/Code/lb/CMakeLists.txt @@ -16,4 +16,5 @@ add_library(hemelb_lb OBJECT kernels/CassonRheologyModel.cc kernels/TruncatedPowerLawRheologyModel.cc MacroscopicPropertyCache.cc SimulationState.cc StabilityTester.cc InitialCondition.cc + lb.cc ) diff --git a/Code/lb/EntropyTester.h b/Code/lb/EntropyTester.h index f5e9217b2..9b77ddf04 100644 --- a/Code/lb/EntropyTester.h +++ b/Code/lb/EntropyTester.h @@ -11,219 +11,174 @@ #include "lb/HFunction.h" #include "log/Logger.h" -namespace hemelb +namespace hemelb::lb { - namespace lb - { - template + template class EntropyTester : public net::PhasedBroadcastRegular { - public: + private: + enum HTHEOREM + { + OBEYED, + DISOBEYED + }; + + /** + * Slightly arbitrary spread factor for the tree. + */ + static constexpr unsigned SPREADFACTOR = 10; + + geometry::Domain const* mLatDat; + + /** + * Stability value of this node and its children to propagate upwards. + */ + int mUpwardsValue = OBEYED; + + static constexpr auto ObeyedArray() { + std::array x; + std::fill(x.begin(), x.end(), OBEYED); + return x; + } + + /** + * Array for storing the passed-up stability values from child nodes. + */ + std::array mChildrensValues = ObeyedArray(); + + std::array mCollisionTypesTested; + std::unique_ptr mHPreCollision; + + public: EntropyTester(int* collisionTypes, unsigned int typesTested, const geometry::Domain * iLatDat, net::Net* net, SimulationState* simState) : - net::PhasedBroadcastRegular(net, simState, SPREADFACTOR), + net::PhasedBroadcastRegular(net, simState, SPREADFACTOR), mLatDat(iLatDat) { - for (unsigned int i = 0; i < COLLISION_TYPES; i++) - { - mCollisionTypesTested[i] = false; - } - for (unsigned int i = 0; i < typesTested; i++) - { - mCollisionTypesTested[collisionTypes[i]] = true; - } - - mHPreCollision = new double[mLatDat->GetLocalFluidSiteCount()]; - - Reset(); - } + std::fill(mCollisionTypesTested.begin(), mCollisionTypesTested.end(), false); + for (unsigned int i = 0; i < typesTested; i++) + { + mCollisionTypesTested[collisionTypes[i]] = true; + } - ~EntropyTester() - { - delete[] mHPreCollision; + mHPreCollision.reset(new double[mLatDat->GetLocalFluidSiteCount()]); + + Reset(); } - void PreReceive() + ~EntropyTester() override = default; + + void PreReceive() override { - site_t offset = 0; - double dHMax = 0.0; - - // The order of arguments in max is important - // If distributions go negative HFunc.eval() will return a NaN. Due to the - // nature of NaN and the structure of max, dH will be assigned as NaN if this is the case - // This is what we want, because the EntropyTester will fail otherwise and abort when - // it is simply sufficient to wait until StabilityTester restarts. - - for (unsigned int collision_type = 0; collision_type < COLLISION_TYPES; collision_type++) - { - if (mCollisionTypesTested[collision_type]) - { - for (site_t i = offset; - i < offset + mLatDat->GetMidDomainCollisionCount(collision_type); i++) - { - const geometry::Site site = mLatDat->GetSite(i); - - HFunction HFunc(site.GetFOld(), nullptr); - dHMax = std::max(dHMax, HFunc.eval() - mHPreCollision[i]); - } + double dHMax = 0.0; + + // The order of arguments in max is important + // If distributions go negative HFunc.eval() will return a NaN. Due to the + // nature of NaN and the structure of max, dH will be assigned as NaN if this is the case + // This is what we want, because the EntropyTester will fail otherwise and abort when + // it is simply sufficient to wait until StabilityTester restarts. + + for (unsigned collision_type = 0; collision_type < COLLISION_TYPES; collision_type++) { + if (mCollisionTypesTested[collision_type]) { + auto [begin, end] = mLatDat->GetMidDomainSiteRange(collision_type); + for (site_t i = begin; i < end; ++i) { + auto const site = mLatDat->GetSite(i); + + HFunction HFunc(site.GetFOld(), nullptr); + dHMax = std::max(dHMax, HFunc.eval() - mHPreCollision[i]); + } + } } - offset += mLatDat->GetMidDomainCollisionCount(collision_type); - } + for (unsigned collision_type = 0; collision_type < COLLISION_TYPES; collision_type++) { + auto [begin, end] = mLatDat->GetDomainEdgeSiteRange(collision_type); + if (mCollisionTypesTested[collision_type]) { + for (site_t i = begin; i < end; ++i) { + auto const site = mLatDat->GetSite(i); - for (unsigned int collision_type = 0; collision_type < COLLISION_TYPES; collision_type++) - { - if (mCollisionTypesTested[collision_type]) - { - for (site_t i = offset; - i < offset + mLatDat->GetDomainEdgeCollisionCount(collision_type); i++) - { - const geometry::Site site = mLatDat->GetSite(i); - - HFunction HFunc(site.GetFOld(), nullptr); - dHMax = std::max(dHMax, HFunc.eval() - mHPreCollision[i]); - } + HFunction HFunc(site.GetFOld(), nullptr); + dHMax = std::max(dHMax, HFunc.eval() - mHPreCollision[i]); + } + } } - offset += mLatDat->GetDomainEdgeCollisionCount(collision_type); - } - - /* - * Ideally dH should never be greater than zero. However, because H is positive definite accuracy as well as - * rounding and truncation errors can make dH greater than zero in certain cases. The tolerance - * is limited by the accuracy of the simulation (including the accuracy to which alpha is calculated) - * The tolerance has to be at least as big as the accuracy to which alpha is calculated - */ - if (dHMax > 1.0E-6) - { - mUpwardsValue = DISOBEYED; - } + // Ideally dH should never be greater than zero. However, because H is positive + // definite accuracy as well as rounding and truncation errors can make dH greater + // than zero in certain cases. The tolerance is limited by the accuracy of the + // simulation (including the accuracy to which alpha is calculated) The tolerance + // has to be at least as big as the accuracy to which alpha is calculated. + if (dHMax > 1.0E-6) + mUpwardsValue = DISOBEYED; } - /** - * Override the reset method in the base class, to reset the stability variables. - */ + // Override the reset method in the base class, to reset the stability variables. void Reset() { - // Re-initialise all values to indicate that the H-theorem is obeyed. - mUpwardsValue = OBEYED; - - for (unsigned int ii = 0; ii < SPREADFACTOR; ii++) - { - mChildrensValues[ii] = OBEYED; - } + // Re-initialise all values to indicate that the H-theorem is obeyed. + mUpwardsValue = OBEYED; + mChildrensValues = ObeyedArray(); } - protected: - /** - * Override the methods from the base class to propagate data from the root, and - * to send data about this node and its childrens' stabilities up towards the root. - */ - void ProgressFromChildren(unsigned long splayNumber) + protected: + // Override the methods from the base class to propagate data from the root, and + // to send data about stability of this node and its children up towards the root. + void ProgressFromChildren(unsigned long splayNumber) override { - ReceiveFromChildren(mChildrensValues, 1); + ReceiveFromChildren(mChildrensValues.data(), 1); } - void ProgressToParent(unsigned long splayNumber) + void ProgressToParent(unsigned long splayNumber) override { - // Store pre-collision values. - site_t offset = 0; - - for (unsigned int collision_type = 0; collision_type < COLLISION_TYPES; collision_type++) - { - if (mCollisionTypesTested[collision_type]) - { - for (site_t i = offset; - i < offset + mLatDat->GetMidDomainCollisionCount(collision_type); i++) - { - const geometry::Site site = mLatDat->GetSite(i); - HFunction HFunc(site.GetFOld(), nullptr); - mHPreCollision[i] = HFunc.eval(); - } + for (unsigned collision_type = 0; collision_type < COLLISION_TYPES; collision_type++) { + if (mCollisionTypesTested[collision_type]) { + auto [begin, end] = mLatDat->GetMidDomainSiteRange(collision_type); + for (site_t i = begin; i < end; ++i) { + auto const site = mLatDat->GetSite(i); + HFunction HFunc(site.GetFOld(), nullptr); + mHPreCollision[i] = HFunc.eval(); + } + } } - offset += mLatDat->GetMidDomainCollisionCount(collision_type); - } - - for (unsigned int collision_type = 0; collision_type < COLLISION_TYPES; collision_type++) - { - if (mCollisionTypesTested[collision_type]) - { - for (site_t i = offset; - i < offset + mLatDat->GetDomainEdgeCollisionCount(collision_type); i++) - { - const geometry::Site site = mLatDat->GetSite(i); - HFunction HFunc(site.GetFOld(), nullptr); - mHPreCollision[i] = HFunc.eval(); - } + for (unsigned collision_type = 0; collision_type < COLLISION_TYPES; collision_type++) { + if (mCollisionTypesTested[collision_type]) { + auto [begin, end] = mLatDat->GetDomainEdgeSiteRange(collision_type); + for (site_t i = begin; i < end; ++i) { + auto const site = mLatDat->GetSite(i); + HFunction HFunc(site.GetFOld(), nullptr); + mHPreCollision[i] = HFunc.eval(); + } + } } - - offset += mLatDat->GetDomainEdgeCollisionCount(collision_type); - } - - SendToParent(&mUpwardsValue, 1); + SendToParent(&mUpwardsValue, 1); } /** * Take the combined stability information (an int, with a value of hemelb::lb::Unstable * if any child node is unstable) and start passing it back down the tree. */ - void TopNodeAction() - { - if (mUpwardsValue == DISOBEYED) - { - log::Logger::Log("H Theorem violated."); - } - } - - /** - * Override the method from the base class to use the data from child nodes. - */ - void PostReceiveFromChildren(unsigned long splayNumber) + void TopNodeAction() override { - // No need to test children's entropy direction if this node already disobeys H-theorem. - if (mUpwardsValue == OBEYED) - { - for (int ii = 0; ii < (int) SPREADFACTOR; ii++) + if (mUpwardsValue == DISOBEYED) { - if (mChildrensValues[ii] == DISOBEYED) - { - mUpwardsValue = DISOBEYED; - break; - } + log::Logger::Log("H Theorem violated."); } - } } - private: - enum HTHEOREM - { - OBEYED, - DISOBEYED - }; - - /** - * Slightly arbitrary spread factor for the tree. - */ - static const unsigned int SPREADFACTOR = 10; - - const geometry::Domain * mLatDat; - - /** - * Stability value of this node and its children to propagate upwards. - */ - int mUpwardsValue; - /** - * Array for storing the passed-up stability values from child nodes. - */ - int mChildrensValues[SPREADFACTOR]; - - bool mCollisionTypesTested[COLLISION_TYPES]; - double* mHPreCollision; + // Override the method from the base class to use the data from child nodes. + void PostReceiveFromChildren(unsigned long splayNumber) override { + // No need to test children's entropy direction if this node already disobeys H-theorem. + if (mUpwardsValue == OBEYED) { + if (std::any_of( + mChildrensValues.begin(), mChildrensValues.end(), + [](int _) { return _ == DISOBEYED; } + )) { + mUpwardsValue = DISOBEYED; + } + } + } }; - - } } #endif /* HEMELB_LB_ENTROPYTESTER_H */ diff --git a/Code/lb/LbmParameters.h b/Code/lb/LbmParameters.h index 3df4cea3c..325bb0f72 100644 --- a/Code/lb/LbmParameters.h +++ b/Code/lb/LbmParameters.h @@ -132,9 +132,6 @@ namespace hemelb::lb // Assume the first site to be used in the kernel is the first site in the core, unless otherwise specified InitParams() = default; - // The number of sites using this kernel instance. - site_t siteCount; - // Each streamer is responsible for updating certain types of sites. These are arranged such they are largely // contiguous in memory (the local contiguous site id). This data structure refers to which of those are handled // by the current streamer. These are given as a collection of contiguous site ids, running from e.g. diff --git a/Code/lb/lb.cc b/Code/lb/lb.cc new file mode 100644 index 000000000..3ee39304e --- /dev/null +++ b/Code/lb/lb.cc @@ -0,0 +1,73 @@ +// This file is part of HemeLB and is Copyright (C) +// the HemeLB team and/or their institutions, as detailed in the +// file AUTHORS. This software is provided under the terms of the +// license in the file LICENSE. + +#include "lb/lb.hpp" + +namespace hemelb::lb { + + LBMBase::LBMBase(LbmParameters iParams, net::Net* net, + geometry::FieldData* latDat, SimulationState* simState, + reporting::Timers &atimings, + geometry::neighbouring::NeighbouringDataManager *neighbouringDataManager) : + mNet(net), mLatDat(latDat), mState(simState), + mParams(std::move(iParams)), + timings(atimings), propertyCache(*simState, latDat->GetDomain()), + neighbouringDataManager(neighbouringDataManager) + { + } + + void LBMBase::Initialise(BoundaryValues* iInletValues, + BoundaryValues* iOutletValues) + { + mInletValues = iInletValues; + mOutletValues = iOutletValues; + + InitCollisions(); + } + + void LBMBase::PrepareBoundaryObjects() + { + // First, iterate through all of the inlet and outlet objects, finding out the minimum density seen in the simulation. + distribn_t minDensity = std::numeric_limits::max(); + + for (unsigned inlet = 0; inlet < mInletValues->GetLocalIoletCount(); ++inlet) + { + minDensity = std::min(minDensity, mInletValues->GetLocalIolet(inlet)->GetDensityMin()); + } + + for (unsigned outlet = 0; outlet < mOutletValues->GetLocalIoletCount(); ++outlet) + { + minDensity = std::min(minDensity, mOutletValues->GetLocalIolet(outlet)->GetDensityMin()); + } + + // Now go through them again, informing them of the minimum density. + for (unsigned inlet = 0; inlet < mInletValues->GetLocalIoletCount(); ++inlet) + { + mInletValues->GetLocalIolet(inlet)->SetMinimumSimulationDensity(minDensity); + } + + for (unsigned outlet = 0; outlet < mOutletValues->GetLocalIoletCount(); ++outlet) + { + mOutletValues->GetLocalIolet(outlet)->SetMinimumSimulationDensity(minDensity); + } + } + + void LBMBase::RequestComms() + { + timings.lb().Start(); + + // Delegate to the lattice data object to post the asynchronous sends and receives + // (via the Net object). + // NOTE that this doesn't actually *perform* the sends and receives, it asks the Net + // to include them in the ISends and IRecvs that happen later. + mLatDat->SendAndReceive(mNet); + + timings.lb().Stop(); + } + + void LBMBase::EndIteration() + { + } +} \ No newline at end of file diff --git a/Code/lb/lb.h b/Code/lb/lb.h index 55c4d8e48..0f6c5d7f2 100644 --- a/Code/lb/lb.h +++ b/Code/lb/lb.h @@ -23,17 +23,72 @@ */ namespace hemelb::lb { + + /// Holds LBM data and behaviour that does not depend on the velocity set/collision etc. + class LBMBase : public net::IteratedAction { + protected: + net::Net* mNet; + geometry::FieldData* mLatDat; + SimulationState* mState; + BoundaryValues *mInletValues = nullptr; + BoundaryValues* mOutletValues = nullptr; + + LbmParameters mParams; + + hemelb::reporting::Timers &timings; + + MacroscopicPropertyCache propertyCache; + + geometry::neighbouring::NeighbouringDataManager *neighbouringDataManager; + + public: + /// OK, with apologies we have a 2-stage construction process. + /// The result of this constructor isn't yet ready to simulate, it must + /// have Initialise(...) called also. + /// + /// Constructor separated due to need to access the partially initialized + /// LBM in order to initialize the arguments to the second construction phase. + LBMBase(LbmParameters params, net::Net* net, + geometry::FieldData* latDat, SimulationState* simState, reporting::Timers &atimings, + geometry::neighbouring::NeighbouringDataManager *neighbouringDataManager); + + /// Second constructor. + void Initialise(BoundaryValues* iInletValues, + BoundaryValues* iOutletValues); + + void RequestComms() override; ///< part of IteratedAction interface. + void EndIteration() override; ///< part of IteratedAction interface. + + inline LbmParameters* GetLbmParams() { + return &mParams; + } + inline MacroscopicPropertyCache& GetPropertyCache() { + return propertyCache; + } + + /// Interface to apply initial conditions + virtual void SetInitialConditions(lb::InitialCondition const& ic_conf, const net::IOCommunicator& ioComms) = 0; + + protected: + virtual void InitCollisions() = 0; + + /** + * Ensure that the BoundaryValues objects have all necessary fields populated. + */ + void PrepareBoundaryObjects(); + }; + /** * Class providing core Lattice Boltzmann functionality. * Implements the IteratedAction interface. */ template> - class LBM : public net::IteratedAction + class LBM : public LBMBase { - public: + public: //! Instantiation type using Traits = TRAITS; - private: + private: using LatticeType = typename Traits::Lattice; // Use the kernel specified through the build system. This will select one of the above classes. using LBKernel = typename Traits::Kernel; @@ -49,97 +104,45 @@ namespace hemelb::lb using tInletWallCollision = typename Traits::WallInletBoundary; using tOutletWallCollision = typename Traits::WallOutletBoundary; - public: - /** - * Constructor, stage 1. - * Object so initialized is not ready for simulation. - * Must have Initialise(...) called also. Constructor separated due to need to access - * the partially initialized LBM in order to initialize the arguments to the second construction phase. - */ - LBM(LbmParameters params, net::Net* net, - geometry::FieldData* latDat, SimulationState* simState, reporting::Timers &atimings, - geometry::neighbouring::NeighbouringDataManager *neighbouringDataManager); + public: + using LBMBase::LBMBase; ~LBM() override = default; - void RequestComms() override; ///< part of IteratedAction interface. + void PreSend() override; ///< part of IteratedAction interface. void PreReceive() override; ///< part of IteratedAction interface. void PostReceive() override; ///< part of IteratedAction interface. - void EndIteration() override; ///< part of IteratedAction interface. - - [[nodiscard]] site_t TotalFluidSiteCount() const; - void SetTotalFluidSiteCount(site_t); - /** - * Second constructor. - * - */ - void Initialise(BoundaryValues* iInletValues, - BoundaryValues* iOutletValues); - void SetInitialConditions(lb::InitialCondition const& ic_conf, const net::IOCommunicator& ioComms); - - hemelb::lb::LbmParameters *GetLbmParams(); - lb::MacroscopicPropertyCache& GetPropertyCache(); - - private: - - void InitCollisions(); - // The following function pair simplify initialising the site ranges for each collider object. - void InitInitParamsSiteRanges(InitParams& initParams, unsigned& state); - void AdvanceInitParamsSiteRanges(InitParams& initParams, unsigned& state); - /** - * Ensure that the BoundaryValues objects have all necessary fields populated. - */ - void PrepareBoundaryObjects(); + void SetInitialConditions(lb::InitialCondition const& ic_conf, const net::IOCommunicator& ioComms) override; - void ReadParameters(); + private: - void handleIOError(int iError); + void InitCollisions() override; - // Collision objects + /// Collision objects + /// We hold these in a tuple to allow heterogeneous iteration over them. // TODO: these should probably be optional but we can't because some of the collisions aren't copyable - std::unique_ptr mMidFluidCollision; - std::unique_ptr mWallCollision; - std::unique_ptr mInletCollision; - std::unique_ptr mOutletCollision; - std::unique_ptr mInletWallCollision; - std::unique_ptr mOutletWallCollision; - - void StreamAndCollide(streamer auto& s, const site_t iFirstIndex, - const site_t iSiteCount) + using CollisionTuple = std::tuple< + std::unique_ptr, + std::unique_ptr, + std::unique_ptr, + std::unique_ptr, + std::unique_ptr, + std::unique_ptr + >; + CollisionTuple mCollisions; + + void StreamAndCollide(streamer auto& s, const site_t beginIndex, const site_t endIndex) { - s.StreamAndCollide(iFirstIndex, - iSiteCount, - &mParams, - *mLatDat, - propertyCache); + s.StreamAndCollide(beginIndex, endIndex, &mParams, *mLatDat, propertyCache); } - void PostStep(streamer auto& s, const site_t iFirstIndex, const site_t iSiteCount) + void PostStep(streamer auto& s, const site_t beginIndex, const site_t endIndex) { - s.PostStep(iFirstIndex, - iSiteCount, - &mParams, - *mLatDat, - propertyCache); - + s.PostStep(beginIndex, endIndex, &mParams, *mLatDat, propertyCache); } - net::Net* mNet; - geometry::FieldData* mLatDat; - SimulationState* mState; - BoundaryValues *mInletValues = nullptr, - *mOutletValues = nullptr; - - LbmParameters mParams; - - hemelb::reporting::Timers &timings; - - MacroscopicPropertyCache propertyCache; - - geometry::neighbouring::NeighbouringDataManager *neighbouringDataManager; }; - } #endif // HEMELB_LB_LB_H diff --git a/Code/lb/lb.hpp b/Code/lb/lb.hpp index 0a945f184..5a93ef5bc 100644 --- a/Code/lb/lb.hpp +++ b/Code/lb/lb.hpp @@ -7,149 +7,75 @@ #define HEMELB_LB_LB_HPP #include "lb/lb.h" +#include #include "lb/InitialCondition.h" #include "lb/InitialCondition.hpp" namespace hemelb::lb { - template - hemelb::lb::LbmParameters* LBM::GetLbmParams() - { - return &mParams; + namespace detail { + // Helper for below using the standard trick of deducing the + // indices from an index_sequence. + template + void tuple_enumerated_foreach_impl(F &&f, std::tuple& t, std::index_sequence is) { + (std::invoke(f, Is, std::get(t)), ...); + } } - template - lb::MacroscopicPropertyCache& LBM::GetPropertyCache() - { - return propertyCache; - } - - template - LBM::LBM(LbmParameters iParams, net::Net* net, - geometry::FieldData* latDat, SimulationState* simState, - reporting::Timers &atimings, - geometry::neighbouring::NeighbouringDataManager *neighbouringDataManager) : - mNet(net), mLatDat(latDat), mState(simState), - mParams(std::move(iParams)), - timings(atimings), propertyCache(*simState, latDat->GetDomain()), - neighbouringDataManager(neighbouringDataManager) - { - } - - template - void LBM::InitInitParamsSiteRanges(InitParams& initParams, unsigned& state) - { - auto& dom = mLatDat->GetDomain(); - initParams.siteRanges.resize(2); - - initParams.siteRanges[0].first = 0; - initParams.siteRanges[1].first = dom.GetMidDomainSiteCount(); - state = 0; - initParams.siteRanges[0].second = initParams.siteRanges[0].first - + dom.GetMidDomainCollisionCount(state); - initParams.siteRanges[1].second = initParams.siteRanges[1].first - + dom.GetDomainEdgeCollisionCount(state); - - initParams.siteCount = dom.GetMidDomainCollisionCount(state) - + dom.GetDomainEdgeCollisionCount(state); - } - - template - void LBM::AdvanceInitParamsSiteRanges(InitParams& initParams, unsigned& state) - { - auto& dom = mLatDat->GetDomain(); - initParams.siteRanges[0].first += dom.GetMidDomainCollisionCount(state); - initParams.siteRanges[1].first += dom.GetDomainEdgeCollisionCount(state); - ++state; - initParams.siteRanges[0].second = initParams.siteRanges[0].first - + dom.GetMidDomainCollisionCount(state); - initParams.siteRanges[1].second = initParams.siteRanges[1].first - + dom.GetDomainEdgeCollisionCount(state); - - initParams.siteCount = dom.GetMidDomainCollisionCount(state) - + dom.GetDomainEdgeCollisionCount(state); + // Invoke the function object with arguments (i, std::get(tuple)), + // where i is a std::size_t and the tuple is the other argument. + // + // Since element types of a tuple can be heterogeneous, the function + // object needs to be similar to a variant's visitor and have an + // overloaded or generic operator(). + // + // TODO: consider a const overload + template + void tuple_enumerated_foreach(std::tuple& t, F&& func) { + using I = std::make_index_sequence; + detail::tuple_enumerated_foreach_impl(std::forward(func), t, I{}); } template void LBM::InitCollisions() { - /** - * Ensure the boundary objects have all info necessary. - */ - PrepareBoundaryObjects(); - - // TODO Note that the convergence checking is not yet implemented in the - // new boundary condition hierarchy system. - // It'd be nice to do this with something like - // MidFluidCollision = new ConvergenceCheckingWrapper(new WhateverMidFluidCollision()); - - auto initParams = InitParams(); - initParams.latDat = &mLatDat->GetDomain(); - initParams.lbmParams = &mParams; - initParams.neighbouringDataManager = neighbouringDataManager; - - unsigned collId; - InitInitParamsSiteRanges(initParams, collId); - mMidFluidCollision = std::make_unique(initParams); - - AdvanceInitParamsSiteRanges(initParams, collId); - mWallCollision = std::make_unique(initParams); - - AdvanceInitParamsSiteRanges(initParams, collId); - initParams.boundaryObject = mInletValues; - mInletCollision = std::make_unique(initParams); - - AdvanceInitParamsSiteRanges(initParams, collId); - initParams.boundaryObject = mOutletValues; - mOutletCollision = std::make_unique(initParams); - - AdvanceInitParamsSiteRanges(initParams, collId); - initParams.boundaryObject = mInletValues; - mInletWallCollision = std::make_unique(initParams); - - AdvanceInitParamsSiteRanges(initParams, collId); - initParams.boundaryObject = mOutletValues; - mOutletWallCollision = std::make_unique(initParams); - } - - template - void LBM::Initialise(BoundaryValues* iInletValues, - BoundaryValues* iOutletValues) - { - mInletValues = iInletValues; - mOutletValues = iOutletValues; - - InitCollisions(); + // Ensure the boundary objects have all info necessary. + PrepareBoundaryObjects(); + + // TODO Note that the convergence checking is not yet implemented in the + // new boundary condition hierarchy system. + // It'd be nice to do this with something like + // MidFluidCollision = new ConvergenceCheckingWrapper(new WhateverMidFluidCollision()); + + auto& dom = mLatDat->GetDomain(); + auto initParams = InitParams(); + initParams.latDat = &dom; + initParams.lbmParams = &mParams; + initParams.neighbouringDataManager = neighbouringDataManager; + initParams.siteRanges.resize(2); + + auto get_boundary_vals = [this](std::size_t collision_id) -> BoundaryValues* { + if (collision_id < 2) + return nullptr; + if (collision_id % 2) + return mOutletValues; + else + return mInletValues; + }; + + tuple_enumerated_foreach( + mCollisions, + [&] (std::size_t i, std::unique_ptr& collision) { + initParams.siteRanges[0] = dom.GetMidDomainSiteRange(i); + initParams.siteRanges[1] = dom.GetDomainEdgeSiteRange(i); + initParams.boundaryObject = get_boundary_vals(i); + collision = std::make_unique(initParams); + } + ); } - template - void LBM::PrepareBoundaryObjects() - { - // First, iterate through all of the inlet and outlet objects, finding out the minimum density seen in the simulation. - distribn_t minDensity = std::numeric_limits::max(); - - for (unsigned inlet = 0; inlet < mInletValues->GetLocalIoletCount(); ++inlet) - { - minDensity = std::min(minDensity, mInletValues->GetLocalIolet(inlet)->GetDensityMin()); - } - - for (unsigned outlet = 0; outlet < mOutletValues->GetLocalIoletCount(); ++outlet) - { - minDensity = std::min(minDensity, mOutletValues->GetLocalIolet(outlet)->GetDensityMin()); - } - // Now go through them again, informing them of the minimum density. - for (unsigned inlet = 0; inlet < mInletValues->GetLocalIoletCount(); ++inlet) - { - mInletValues->GetLocalIolet(inlet)->SetMinimumSimulationDensity(minDensity); - } - - for (unsigned outlet = 0; outlet < mOutletValues->GetLocalIoletCount(); ++outlet) - { - mOutletValues->GetLocalIolet(outlet)->SetMinimumSimulationDensity(minDensity); - } - } template void LBM::SetInitialConditions(lb::InitialCondition const& icond, const net::IOCommunicator& ioComms) @@ -158,159 +84,84 @@ namespace hemelb::lb icond.SetTime(mState); } - template - void LBM::RequestComms() - { - timings.lb().Start(); - - // Delegate to the lattice data object to post the asynchronous sends and receives - // (via the Net object). - // NOTE that this doesn't actually *perform* the sends and receives, it asks the Net - // to include them in the ISends and IRecvs that happen later. - mLatDat->SendAndReceive(mNet); - - timings.lb().Stop(); - } - template void LBM::PreSend() { - timings.lb().Start(); - timings.lb_calc().Start(); - - /** - * In the PreSend phase, we do LB on all the sites that need to have results sent to - * neighbouring ranks ('domainEdge' sites). In site id terms, this means we start at the - * end of the sites whose neighbours all lie on this rank ('midDomain'), then progress - * through the sites of each type in turn. - */ - auto& dom = mLatDat->GetDomain(); - site_t offset = dom.GetMidDomainSiteCount(); - - log::Logger::Log("LBM - PreSend - StreamAndCollide"); - StreamAndCollide(*mMidFluidCollision, offset, dom.GetDomainEdgeCollisionCount(0)); - offset += dom.GetDomainEdgeCollisionCount(0); - - StreamAndCollide(*mWallCollision, offset, dom.GetDomainEdgeCollisionCount(1)); - offset += dom.GetDomainEdgeCollisionCount(1); - - mInletValues->FinishReceive(); - StreamAndCollide(*mInletCollision, offset, dom.GetDomainEdgeCollisionCount(2)); - offset += dom.GetDomainEdgeCollisionCount(2); - - mOutletValues->FinishReceive(); - StreamAndCollide(*mOutletCollision, offset, dom.GetDomainEdgeCollisionCount(3)); - offset += dom.GetDomainEdgeCollisionCount(3); - - StreamAndCollide(*mInletWallCollision, offset, dom.GetDomainEdgeCollisionCount(4)); - offset += dom.GetDomainEdgeCollisionCount(4); - - StreamAndCollide(*mOutletWallCollision, offset, dom.GetDomainEdgeCollisionCount(5)); - - timings.lb_calc().Stop(); - timings.lb().Stop(); + timings.lb().Start(); + timings.lb_calc().Start(); + + // In the PreSend phase, we do LB on all the sites that need to have results sent to + // neighbouring ranks ('domainEdge' sites). In site id terms, this means we start at the + // end of the sites whose neighbours all lie on this rank ('midDomain'), then progress + // through the sites of each type in turn. + log::Logger::Log("LBM - PreSend - StreamAndCollide"); + tuple_enumerated_foreach( + mCollisions, + [this] (std::size_t i, auto& collision) { + auto [begin, end] = mLatDat->GetDomain().GetDomainEdgeSiteRange(i); + StreamAndCollide(*collision, begin, end); + } + ); + + timings.lb_calc().Stop(); + timings.lb().Stop(); } template void LBM::PreReceive() { - timings.lb().Start(); - timings.lb_calc().Start(); - - /** - * In the PreReceive phase, we perform LB for all the sites whose neighbours lie on this - * rank ('midDomain' rather than 'domainEdge' sites). Ideally this phase is the longest bit (maximising time for the asynchronous sends - * and receives to complete). - * - * In site id terms, this means starting at the first site and progressing through the - * midDomain sites, one type at a time. - */ - auto& dom = mLatDat->GetDomain(); - site_t offset = 0; - - log::Logger::Log("LBM - PreReceive - StreamAndCollide"); - StreamAndCollide(*mMidFluidCollision, offset, dom.GetMidDomainCollisionCount(0)); - offset += dom.GetMidDomainCollisionCount(0); - - StreamAndCollide(*mWallCollision, offset, dom.GetMidDomainCollisionCount(1)); - offset += dom.GetMidDomainCollisionCount(1); - - StreamAndCollide(*mInletCollision, offset, dom.GetMidDomainCollisionCount(2)); - offset += dom.GetMidDomainCollisionCount(2); - - StreamAndCollide(*mOutletCollision, offset, dom.GetMidDomainCollisionCount(3)); - offset += dom.GetMidDomainCollisionCount(3); - - StreamAndCollide(*mInletWallCollision, offset, dom.GetMidDomainCollisionCount(4)); - offset += dom.GetMidDomainCollisionCount(4); - - StreamAndCollide(*mOutletWallCollision, offset, dom.GetMidDomainCollisionCount(5)); - - timings.lb_calc().Stop(); - timings.lb().Stop(); + timings.lb().Start(); + timings.lb_calc().Start(); + + // In the PreReceive phase, we perform LB for all the sites whose neighbours lie on this + // rank ('midDomain' rather than 'domainEdge' sites). Ideally this phase is the longest bit (maximising time for the asynchronous sends + // and receives to complete). + // + // In site id terms, this means starting at the first site and progressing through the + // midDomain sites, one type at a time. + log::Logger::Log("LBM - PreReceive - StreamAndCollide"); + tuple_enumerated_foreach( + mCollisions, + [this] (std::size_t i, auto& collision) { + auto [begin, end] = mLatDat->GetDomain().GetMidDomainSiteRange(i); + StreamAndCollide(*collision, begin, end); + } + ); + timings.lb_calc().Stop(); + timings.lb().Stop(); } template void LBM::PostReceive() { - timings.lb().Start(); - - // Copy the distribution functions received from the neighbouring - // processors into the destination buffer "f_new". - // This is done here, after receiving the sent distributions from neighbours. - mLatDat->CopyReceived(); - - auto& dom = mLatDat->GetDomain(); - // Do any cleanup steps necessary on boundary nodes - site_t offset = dom.GetMidDomainSiteCount(); - - timings.lb_calc().Start(); - - log::Logger::Log("LBM - PostReceive - StreamAndCollide"); - //TODO yup, this is horrible. If you read this, please improve the following code. - PostStep(*mMidFluidCollision, offset, dom.GetDomainEdgeCollisionCount(0)); - offset += dom.GetDomainEdgeCollisionCount(0); - - PostStep(*mWallCollision, offset, dom.GetDomainEdgeCollisionCount(1)); - offset += dom.GetDomainEdgeCollisionCount(1); - - PostStep(*mInletCollision, offset, dom.GetDomainEdgeCollisionCount(2)); - offset += dom.GetDomainEdgeCollisionCount(2); - - PostStep(*mOutletCollision, offset, dom.GetDomainEdgeCollisionCount(3)); - offset += dom.GetDomainEdgeCollisionCount(3); - - PostStep(*mInletWallCollision, offset, dom.GetDomainEdgeCollisionCount(4)); - offset += dom.GetDomainEdgeCollisionCount(4); - - PostStep(*mOutletWallCollision, offset, dom.GetDomainEdgeCollisionCount(5)); - - offset = 0; - - PostStep(*mMidFluidCollision, offset, dom.GetMidDomainCollisionCount(0)); - offset += dom.GetMidDomainCollisionCount(0); - - PostStep(*mWallCollision, offset, dom.GetMidDomainCollisionCount(1)); - offset += dom.GetMidDomainCollisionCount(1); - - PostStep(*mInletCollision, offset, dom.GetMidDomainCollisionCount(2)); - offset += dom.GetMidDomainCollisionCount(2); - - PostStep(*mOutletCollision, offset, dom.GetMidDomainCollisionCount(3)); - offset += dom.GetMidDomainCollisionCount(3); - - PostStep(*mInletWallCollision, offset, dom.GetMidDomainCollisionCount(4)); - offset += dom.GetMidDomainCollisionCount(4); - - PostStep(*mOutletWallCollision, offset, dom.GetMidDomainCollisionCount(5)); - - timings.lb_calc().Stop(); - timings.lb().Stop(); - } - - template - void LBM::EndIteration() - { + timings.lb().Start(); + + // Copy the distribution functions received from the neighbouring + // processors into the destination buffer "f_new". + // This is done here, after receiving the sent distributions from neighbours. + mLatDat->CopyReceived(); + + // Do any cleanup steps necessary on boundary nodes + timings.lb_calc().Start(); + + log::Logger::Log("LBM - PostReceive - StreamAndCollide"); + tuple_enumerated_foreach( + mCollisions, + [this] (std::size_t i, auto& collision) { + auto [begin, end] = mLatDat->GetDomain().GetDomainEdgeSiteRange(i); + PostStep(*collision, begin, end); + } + ); + tuple_enumerated_foreach( + mCollisions, + [this] (std::size_t i, auto& collision) { + auto [begin, end] = mLatDat->GetDomain().GetMidDomainSiteRange(i); + PostStep(*collision, begin, end); + } + ); + + timings.lb_calc().Stop(); + timings.lb().Stop(); } } diff --git a/Code/lb/streamers/BulkStreamer.h b/Code/lb/streamers/BulkStreamer.h index 47e3a4332..68cb4953d 100644 --- a/Code/lb/streamers/BulkStreamer.h +++ b/Code/lb/streamers/BulkStreamer.h @@ -64,12 +64,12 @@ namespace hemelb::lb { } - void StreamAndCollide(const site_t firstIndex, const site_t siteCount, + void StreamAndCollide(const site_t beginIndex, const site_t endIndex, const LbmParameters* lbmParams, geometry::FieldData& latDat, lb::MacroscopicPropertyCache& propertyCache) { - for (site_t siteIndex = firstIndex; siteIndex < (firstIndex + siteCount); siteIndex++) + for (site_t siteIndex = beginIndex; siteIndex < endIndex; ++siteIndex) { geometry::Site site = latDat.GetSite(siteIndex); VarsType hydroVars(site); @@ -90,7 +90,7 @@ namespace hemelb::lb } } - void PostStep(const site_t iFirstIndex, const site_t iSiteCount, + void PostStep(const site_t beginIndex, const site_t iSiteCount, const LbmParameters* iLbmParams, geometry::FieldData& bLatDat, lb::MacroscopicPropertyCache& propertyCache) { diff --git a/Code/lb/streamers/JunkYang.h b/Code/lb/streamers/JunkYang.h index 97549ee22..44fb194e2 100644 --- a/Code/lb/streamers/JunkYang.h +++ b/Code/lb/streamers/JunkYang.h @@ -76,12 +76,12 @@ namespace hemelb::lb } } - void StreamAndCollide(const site_t firstIndex, const site_t siteCount, + void StreamAndCollide(const site_t firstIndex, const site_t lastIndex, const LbmParameters* lbmParams, geometry::FieldData& latticeData, lb::MacroscopicPropertyCache& propertyCache) { - for (site_t siteIndex = firstIndex; siteIndex < (firstIndex + siteCount); siteIndex++) + for (site_t siteIndex = firstIndex; siteIndex < lastIndex; siteIndex++) { HASSERT(latticeData.GetSite(siteIndex).IsWall()); HASSERT(lMatrices.find(siteIndex) != lMatrices.end()); @@ -153,11 +153,11 @@ namespace hemelb::lb } - void PostStep(const site_t firstIndex, const site_t siteCount, + void PostStep(const site_t firstIndex, const site_t lastIndex, const LbmParameters* lbmParams, geometry::FieldData& latticeData, lb::MacroscopicPropertyCache& propertyCache) { - for (site_t siteIndex = firstIndex; siteIndex < (firstIndex + siteCount); siteIndex++) + for (site_t siteIndex = firstIndex; siteIndex < lastIndex; siteIndex++) { HASSERT(latticeData.GetSite(siteIndex).IsWall()); HASSERT(lMatrices.find(siteIndex) != lMatrices.end()); diff --git a/Code/lb/streamers/StreamerTypeFactory.h b/Code/lb/streamers/StreamerTypeFactory.h index 30a7517f2..ba2706ca3 100644 --- a/Code/lb/streamers/StreamerTypeFactory.h +++ b/Code/lb/streamers/StreamerTypeFactory.h @@ -45,12 +45,12 @@ namespace hemelb::lb { } - void StreamAndCollide(const site_t firstIndex, const site_t siteCount, + void StreamAndCollide(const site_t beginIndex, const site_t endIndex, const LbmParameters* lbmParams, geometry::FieldData& latDat, lb::MacroscopicPropertyCache& propertyCache) { - for (site_t siteIndex = firstIndex; siteIndex < (firstIndex + siteCount); siteIndex++) + for (site_t siteIndex = beginIndex; siteIndex < endIndex; ++siteIndex) { geometry::Site site = latDat.GetSite(siteIndex); VarsType hydroVars(site); @@ -85,11 +85,11 @@ namespace hemelb::lb } } - void PostStep(const site_t firstIndex, const site_t siteCount, + void PostStep(const site_t beginIndex, const site_t endIndex, const LbmParameters* lbmParams, geometry::FieldData& latticeData, lb::MacroscopicPropertyCache& propertyCache) { - for (site_t siteIndex = firstIndex; siteIndex < (firstIndex + siteCount); siteIndex++) + for (site_t siteIndex = beginIndex; siteIndex < endIndex; ++siteIndex) { geometry::Site site = latticeData.GetSite(siteIndex); for (unsigned int direction = 0; direction < LatticeType::NUMVECTORS; direction++) diff --git a/Code/lb/streamers/VirtualSiteIolet.h b/Code/lb/streamers/VirtualSiteIolet.h index c62115c52..d07ae7cfd 100644 --- a/Code/lb/streamers/VirtualSiteIolet.h +++ b/Code/lb/streamers/VirtualSiteIolet.h @@ -130,13 +130,13 @@ namespace hemelb::lb * links will be done in the post-step as we must ensure that all * the data is available to construct virtual sites. */ - void StreamAndCollide(const site_t firstIndex, const site_t siteCount, + void StreamAndCollide(const site_t firstIndex, const site_t lastIndex, const LbmParameters* lbmParams, geometry::FieldData& latDat, lb::MacroscopicPropertyCache& propertyCache) { auto&& dom = latDat.GetDomain(); - for (site_t siteIndex = firstIndex; siteIndex < (firstIndex + siteCount); siteIndex++) + for (site_t siteIndex = firstIndex; siteIndex < lastIndex; siteIndex++) { auto&& site = latDat.GetSite(siteIndex); VarsType hydroVars(site); @@ -184,14 +184,14 @@ namespace hemelb::lb } } - void PostStep(const site_t firstIndex, const site_t siteCount, + void PostStep(const site_t firstIndex, const site_t lastIndex, const LbmParameters* lbmParams, geometry::FieldData* latDat, lb::MacroscopicPropertyCache& propertyCache) { const LatticeTimeStep t = bValues->GetTimeStep(); const typename VSiteByLocalIdxMultiMap::iterator beginVSites = vsByLocalIdx.lower_bound(firstIndex), endVSites = - vsByLocalIdx.lower_bound(firstIndex + siteCount); + vsByLocalIdx.lower_bound(lastIndex); for (typename VSiteByLocalIdxMultiMap::iterator vSiteIt = beginVSites; vSiteIt != endVSites; ++vSiteIt) diff --git a/Code/multiscale/MultiscaleSimulationMaster.h b/Code/multiscale/MultiscaleSimulationMaster.h index 04c3413c2..aef65761b 100644 --- a/Code/multiscale/MultiscaleSimulationMaster.h +++ b/Code/multiscale/MultiscaleSimulationMaster.h @@ -21,7 +21,7 @@ namespace hemelb::multiscale template class MultiscaleSimulationMaster : public SimulationMaster<> { - public: + public: MultiscaleSimulationMaster(configuration::CommandLine &options, const net::IOCommunicator& ioComm, Intercommunicator & aintercomms) : @@ -87,39 +87,26 @@ namespace hemelb::multiscale if (velocity == true) { /* Do not include the non-iolet adjacent sites (resp. MidFluid and Wall-adjacent). */ - long long int offset = domainData->GetMidDomainCollisionCount(0) - + domainData->GetMidDomainCollisionCount(1); /* Do include iolet adjacent sites (inlet) */ - long long int ioletsSiteCount = domainData->GetMidDomainCollisionCount(2); invertedInletBoundaryList = PopulateInvertedBoundaryList(domainData.get(), invertedInletBoundaryList, - offset, - ioletsSiteCount); + domainData->GetMidDomainSiteRange(2)); - offset += domainData->GetMidDomainCollisionCount(2); /* Do include iolet adjacent sites (outlet) */ - ioletsSiteCount = domainData->GetMidDomainCollisionCount(3); invertedOutletBoundaryList = PopulateInvertedBoundaryList(domainData.get(), invertedOutletBoundaryList, - offset, - ioletsSiteCount); + domainData->GetMidDomainSiteRange(3)); - offset += domainData->GetMidDomainCollisionCount(3); /* Do include iolet adjacent sites (inlet-wall) */ - ioletsSiteCount = domainData->GetMidDomainCollisionCount(4); invertedInletBoundaryList = PopulateInvertedBoundaryList(domainData.get(), invertedInletBoundaryList, - offset, - ioletsSiteCount); + domainData->GetMidDomainSiteRange(4)); - offset += domainData->GetMidDomainCollisionCount(4); /* Do include iolet adjacent sites (outlet-wall) */ - ioletsSiteCount = domainData->GetMidDomainCollisionCount(5); invertedOutletBoundaryList = PopulateInvertedBoundaryList(domainData.get(), invertedOutletBoundaryList, - offset, - ioletsSiteCount); + domainData->GetMidDomainSiteRange(5)); log::Logger::Log("Populated inlets (numinlets/sizeinlet0): %i/%i", invertedInletBoundaryList.size(), @@ -282,7 +269,7 @@ namespace hemelb::multiscale Intercommunicator &intercomms; typename Intercommunicator::IntercommunicandTypeT multiscaleIoletType; - private: + private: /* Loops over iolets to set the need for communications. */ void SetCommsRequired(lb::BoundaryValues* ioletValues, bool b) @@ -309,19 +296,18 @@ namespace hemelb::multiscale //out: iBL //in: domain_type, [Site Object]->SiteData->GetBoundaryID, //MPI_Gatherv if data is only this process. - std::vector > PopulateInvertedBoundaryList( + std::vector > + PopulateInvertedBoundaryList( geometry::Domain* latticeData, - std::vector > invertedBoundaryList, int offset, int ioletsSiteCount) - { - for (int i = 0; i < ioletsSiteCount; i++) - { - /* 1. Obtain Boundary ID number. */ - auto&& s = latticeData->GetSite(offset + i).GetSiteData(); - int boundaryID = s.GetIoletId(); - - /* 2. Grow the list to an appropriate size if needed. */ - while ( ((int) invertedBoundaryList.size()) <= boundaryID) - { + std::vector > invertedBoundaryList, std::pair i_range + ) { + for (auto i = i_range.first; i < i_range.second; ++i) { + /* 1. Obtain Boundary ID number. */ + auto&& s = latticeData->GetSite(i).GetSiteData(); + int boundaryID = s.GetIoletId(); + + /* 2. Grow the list to an appropriate size if needed. */ + while (std::ssize(invertedBoundaryList) <= boundaryID) { log::Logger::Log("WARNING: Growing the invertedBoundaryList, because we created in wrongly in MultiscaleSimulation Master."); std::vector a(0); invertedBoundaryList.push_back(a); @@ -330,11 +316,14 @@ namespace hemelb::multiscale log::Logger::Log("ERROR: invertedBoundaryList is growing to ridiculous proportions due to a faulty boundaryID."); exit(-1); } - } - - /* 3. Insert this site in the inverted Boundary List. */ - invertedBoundaryList[boundaryID].push_back(latticeData->GetGlobalNoncontiguousSiteIdFromGlobalCoords(latticeData->GetSite(offset - + i).GetGlobalSiteCoords())); + } + + /* 3. Insert this site in the inverted Boundary List. */ + invertedBoundaryList[boundaryID].push_back( + latticeData->GetGlobalNoncontiguousSiteIdFromGlobalCoords( + latticeData->GetSite(i).GetGlobalSiteCoords() + ) + ); } return invertedBoundaryList; } diff --git a/Code/tests/helpers/FourCubeBasedTestFixture.cc b/Code/tests/helpers/FourCubeBasedTestFixture.cc index 6cdcca1c7..12f6cbc52 100644 --- a/Code/tests/helpers/FourCubeBasedTestFixture.cc +++ b/Code/tests/helpers/FourCubeBasedTestFixture.cc @@ -27,7 +27,6 @@ namespace hemelb::tests::helpers unitConverter = simBuilder->GetUnitConverter(); initParams.latDat = &latDat->GetDomain(); - initParams.siteCount = initParams.latDat->GetLocalFluidSiteCount(); initParams.lbmParams = &lbmParams; numSites = initParams.latDat->GetLocalFluidSiteCount(); } diff --git a/Code/tests/lb/GuoForcingTests.cc b/Code/tests/lb/GuoForcingTests.cc index 7d1fdbdc0..56996018a 100644 --- a/Code/tests/lb/GuoForcingTests.cc +++ b/Code/tests/lb/GuoForcingTests.cc @@ -313,6 +313,7 @@ namespace hemelb::tests // Same for forces LatticeVector const position(2, 2, 2); auto&& site = latDat->GetSite(position); + auto i_site = site.GetIndex(); helpers::allZeroButOne(*latDat, position); // Get collided distributions at this site @@ -325,7 +326,7 @@ namespace hemelb::tests // Stream that site using lb::BulkStreamer; BulkStreamer streamer(initParams); - streamer.StreamAndCollide(site.GetIndex(), 1, &lbmParams, *latDat, *propertyCache); + streamer.StreamAndCollide(i_site, i_site + 1, &lbmParams, *latDat, *propertyCache); // Now check streaming worked correctly for (size_t i(0); i < LatticeType::NUMVECTORS; ++i) { @@ -342,6 +343,7 @@ namespace hemelb::tests // Same for forces LatticeVector const position(1, 1, 1); auto&& site = latDat->GetSite(position); + auto i_site = site.GetIndex(); helpers::allZeroButOne(*latDat, position); // Get collided distributions at this site @@ -357,10 +359,10 @@ namespace hemelb::tests lb::NullLink >; SBB streamer(initParams); - streamer.StreamAndCollide(site.GetIndex(), 1, &lbmParams, *latDat, *propertyCache); + streamer.StreamAndCollide(i_site, i_site + 1, &lbmParams, *latDat, *propertyCache); distribn_t const * const actual = helpers::GetFNew(*latDat, position); - bool paranoia(false); + bool paranoia = false; for (size_t i(0); i < LatticeType::NUMVECTORS; ++i) { if (not site.HasWall(i)) continue; diff --git a/Code/tests/lb/StreamerTests.cc b/Code/tests/lb/StreamerTests.cc index d8bb504b6..f7a5b9520 100644 --- a/Code/tests/lb/StreamerTests.cc +++ b/Code/tests/lb/StreamerTests.cc @@ -220,117 +220,101 @@ namespace hemelb::tests } } - SECTION("TestSimpleBounceBack") { - using SBB = StreamerTypeFactory, NullLink>; - // Initialise fOld in the lattice data. We choose values so - // that each site has an anisotropic distribution function, - // and that each site's function is distinguishable. - LbTestsHelper::InitialiseAnisotropicTestData(*latDat); - - site_t firstWallSite = dom->GetMidDomainCollisionCount(0); - site_t wallSitesCount = dom->GetMidDomainCollisionCount(1); - - // Check that the lattice has the expected number of sites - // labeled as pure wall (otherwise this test is void) - REQUIRE(site_t(24) == dom->GetMidDomainCollisionCount(1)); - - site_t offset = 0; - - // Mid-Fluid sites use simple collide and stream - BulkStreamer simpleCollideAndStream(initParams); - - simpleCollideAndStream.StreamAndCollide(offset, - dom->GetMidDomainCollisionCount(0), - &lbmParams, - *latDat, - *propertyCache); - offset += dom->GetMidDomainCollisionCount(0); - - // Wall sites use simple bounce back - SBB simpleBounceBack(initParams); - - simpleBounceBack.StreamAndCollide(offset, - dom->GetMidDomainCollisionCount(1), - &lbmParams, - *latDat, - *propertyCache); - offset += dom->GetMidDomainCollisionCount(1); - - // Consider inlet/outlets and their walls as mid-fluid sites - simpleCollideAndStream.StreamAndCollide(offset, - dom->GetLocalFluidSiteCount() - - offset, - &lbmParams, - *latDat, - *propertyCache); - offset += dom->GetLocalFluidSiteCount() - offset; - - // Sanity check - REQUIRE(offset == dom->GetLocalFluidSiteCount()); - - // Loop over the wall sites and check whether they got - // properly streamed on or bounced back depending on where - // they sit relative to the wall. We ignore mid-Fluid sites - // since StreamAndCollide was tested before. - for (site_t wallSiteLocalIndex = 0; wallSiteLocalIndex < wallSitesCount; wallSiteLocalIndex++) { - site_t streamedToSite = firstWallSite + wallSiteLocalIndex; - const auto streamedSite = latDat->GetSite(streamedToSite); - distribn_t* streamedToFNew = latDat->GetFNew(NUMVECTORS * streamedToSite); - - for (unsigned int streamedDirection = 0; streamedDirection - < NUMVECTORS; ++streamedDirection) { - unsigned oppDirection = LATTICE::INVERSEDIRECTIONS[streamedDirection]; - - // Index of the site streaming to streamedToSite via direction streamedDirection - site_t streamerIndex = streamedSite.GetStreamedIndex (oppDirection); - - // Is streamerIndex a valid index? - if (streamerIndex >= 0 && streamerIndex < (NUMVECTORS * dom->GetLocalFluidSiteCount())) { - // The streamer index is a valid index in the domain, - // therefore stream and collide has happened - site_t streamerSiteId = streamerIndex / NUMVECTORS; - - // Calculate streamerFOld at this site. - distribn_t streamerFOld[NUMVECTORS]; - LbTestsHelper::InitialiseAnisotropicTestData(streamerSiteId, - streamerFOld); - - // Calculate what the value streamed to site streamedToSite should be. - lb::HydroVars streamerHydroVars(streamerFOld); - normalCollision->CalculatePreCollision(streamerHydroVars, streamedSite); - - normalCollision->Collide(&lbmParams, streamerHydroVars); - - // F_new should be equal to the value that was streamed - // from this other site - // in the same direction as we're streaming from. - INFO("BulkStreamer, StreamAndCollide"); - REQUIRE(apprx(streamerHydroVars.GetFPostCollision()[streamedDirection]) == streamedToFNew[streamedDirection]); - } else { - // The streamer index shows that no one has streamed to - // streamedToSite direction streamedDirection, therefore - // bounce back has happened in that site for that - // direction - - // Initialise streamedToSiteFOld with the original data - distribn_t streamerToSiteFOld[NUMVECTORS]; - LbTestsHelper::InitialiseAnisotropicTestData(streamedToSite, - streamerToSiteFOld); - lb::HydroVars hydroVars(streamerToSiteFOld); - normalCollision->CalculatePreCollision(hydroVars, streamedSite); - - // Simulate post-collision using the collision operator. - normalCollision->Collide(&lbmParams, hydroVars); - - // After streaming FNew in a given direction must be f - // post-collision in the opposite direction following - // collision - INFO("Simple bounce-back: site " << streamedToSite << " direction " << streamedDirection); - REQUIRE(streamedToFNew[streamedDirection] == apprx(hydroVars.GetFPostCollision()[oppDirection])); - } - } - } - } + SECTION("TestSimpleBounceBack") { + using SBB = StreamerTypeFactory, NullLink>; + // Initialise fOld in the lattice data. We choose values so + // that each site has an anisotropic distribution function, + // and that each site's function is distinguishable. + LbTestsHelper::InitialiseAnisotropicTestData(*latDat); + + auto [firstWallSite, pastLastWallSite] = dom->GetMidDomainSiteRange(1); + auto wallSitesCount = pastLastWallSite - firstWallSite; + + // Check that the lattice has the expected number of sites + // labeled as pure wall (otherwise this test is void) + REQUIRE(site_t(24) == wallSitesCount); + + // Mid-Fluid sites use simple collide and stream + BulkStreamer simpleCollideAndStream(initParams); + auto [beginBulk, pastLastBulk] = dom->GetMidDomainSiteRange(0); + simpleCollideAndStream.StreamAndCollide( + beginBulk, pastLastBulk, &lbmParams, *latDat, *propertyCache + ); + + // Wall sites use simple bounce back + SBB simpleBounceBack(initParams); + + simpleBounceBack.StreamAndCollide( + firstWallSite, pastLastWallSite, &lbmParams, *latDat, *propertyCache + ); + + // Consider inlet/outlets and their walls as mid-fluid sites + simpleCollideAndStream.StreamAndCollide( + pastLastWallSite,dom->GetLocalFluidSiteCount(), + &lbmParams, *latDat, *propertyCache + ); + + // Loop over the wall sites and check whether they got + // properly streamed on or bounced back depending on where + // they sit relative to the wall. We ignore mid-Fluid sites + // since StreamAndCollide was tested before. + for (site_t i_site = firstWallSite; i_site < pastLastWallSite; ++i_site) { + const auto streamedSite = latDat->GetSite(i_site); + distribn_t* streamedToFNew = latDat->GetFNew(NUMVECTORS * i_site); + + for (unsigned i = 0; i < NUMVECTORS; ++i) { + unsigned oppDirection = LATTICE::INVERSEDIRECTIONS[i]; + + // Index of the site streaming to i_site via direction i + site_t streamerIndex = streamedSite.GetStreamedIndex (oppDirection); + + // Is streamerIndex a valid index? + if (streamerIndex >= 0 && streamerIndex < (NUMVECTORS * dom->GetLocalFluidSiteCount())) { + // The streamer index is a valid index in the domain, + // therefore stream and collide has happened + site_t streamerSiteId = streamerIndex / NUMVECTORS; + + // Calculate streamerFOld at this site. + distribn_t streamerFOld[NUMVECTORS]; + LbTestsHelper::InitialiseAnisotropicTestData(streamerSiteId, + streamerFOld); + + // Calculate what the value streamed to site i_site should be. + lb::HydroVars streamerHydroVars(streamerFOld); + normalCollision->CalculatePreCollision(streamerHydroVars, streamedSite); + + normalCollision->Collide(&lbmParams, streamerHydroVars); + + // F_new should be equal to the value that was streamed + // from this other site + // in the same direction as we're streaming from. + INFO("BulkStreamer, StreamAndCollide"); + REQUIRE(apprx(streamerHydroVars.GetFPostCollision()[i]) == streamedToFNew[i]); + } else { + // The streamer index shows that no one has streamed to + // i_site direction i, therefore + // bounce back has happened in that site for that + // direction + + // Initialise streamedToSiteFOld with the original data + distribn_t streamerToSiteFOld[NUMVECTORS]; + LbTestsHelper::InitialiseAnisotropicTestData(i_site, + streamerToSiteFOld); + lb::HydroVars hydroVars(streamerToSiteFOld); + normalCollision->CalculatePreCollision(hydroVars, streamedSite); + + // Simulate post-collision using the collision operator. + normalCollision->Collide(&lbmParams, hydroVars); + + // After streaming FNew in a given direction must be f + // post-collision in the opposite direction following + // collision + INFO("Simple bounce-back: site " << i_site << " direction " << i); + REQUIRE(streamedToFNew[i] == apprx(hydroVars.GetFPostCollision()[oppDirection])); + } + } + } + } SECTION("GuoZhengShi") { using GZS = StreamerTypeFactory, NullLink>; @@ -496,132 +480,114 @@ namespace hemelb::tests // Junk&Yang should behave like simple bounce back when fluid // sites are 0.5 lattice length units away from the domain // boundary. - SECTION("JunkYangEquivalentToBounceBack") { - // Initialise fOld in the lattice data. We choose values so - // that each site has an anisotropic distribution function, - // and that each site's function is distinguishable. - LbTestsHelper::InitialiseAnisotropicTestData(*latDat); - // Setting all the wall distances to 0.5 will make Junk&Yang - // behave like Simple Bounce Back - LbTestsHelper::SetWallAndIoletDistances(*latDat, 0.5); - - site_t firstWallSite = dom->GetMidDomainCollisionCount(0); - site_t wallSitesCount = dom->GetMidDomainCollisionCount(1) - firstWallSite; - - // Check that the lattice has sites labeled as wall (otherwise - // this test is void) - REQUIRE(wallSitesCount == 16); - - site_t offset = 0; - - // Mid-Fluid sites use simple collide and stream - BulkStreamer simpleCollideAndStream(initParams); - - simpleCollideAndStream.StreamAndCollide(offset, - dom->GetMidDomainCollisionCount(0), - &lbmParams, - *latDat, - *propertyCache); - offset += dom->GetMidDomainCollisionCount(0); - - // Wall sites use Junk and Yang - initParams.siteRanges.push_back(std::pair(offset, - offset - + dom->GetMidDomainCollisionCount(1))); - JunkYangFactory> junkYang(initParams); - - junkYang.StreamAndCollide(offset, - dom->GetMidDomainCollisionCount(1), - &lbmParams, - *latDat, - *propertyCache); - - junkYang.PostStep(offset, - dom->GetMidDomainCollisionCount(1), - &lbmParams, - *latDat, - *propertyCache); - - offset += dom->GetMidDomainCollisionCount(1); - - // Consider inlet/outlets and their walls as mid-fluid sites - simpleCollideAndStream.StreamAndCollide(offset, - dom->GetLocalFluidSiteCount() - - offset, - &lbmParams, - *latDat, - *propertyCache); - offset += dom->GetLocalFluidSiteCount() - offset; - - // Sanity check - REQUIRE(offset == dom->GetLocalFluidSiteCount()); - - // Loop over the wall sites and check whether they got - // properly streamed on or bounced back depending on where - // they sit relative to the wall. We ignore mid-Fluid sites - // since StreamAndCollide was tested before. - for (site_t wallSiteLocalIndex = 0; wallSiteLocalIndex < wallSitesCount; wallSiteLocalIndex++) { - site_t streamedToSite = firstWallSite + wallSiteLocalIndex; - const auto streamedSite = latDat->GetSite(streamedToSite); - distribn_t* streamedToFNew = latDat->GetFNew(NUMVECTORS * streamedToSite); - - for (unsigned int streamedDirection = 0; - streamedDirection < NUMVECTORS; ++streamedDirection) { - unsigned oppDirection = LATTICE::INVERSEDIRECTIONS[streamedDirection]; - - // Index of the site streaming to streamedToSite via - // direction streamedDirection - site_t streamerIndex = streamedSite.GetStreamedIndex (oppDirection); - - // Is streamerIndex a valid index? - if (streamerIndex >= 0 && - streamerIndex < (NUMVECTORS * dom->GetLocalFluidSiteCount())) { - // The streamer index is a valid index in the domain, - // therefore stream and collide has happened - site_t streamerSiteId = streamerIndex / NUMVECTORS; - - // Calculate streamerFOld at this site. - distribn_t streamerFOld[NUMVECTORS]; - LbTestsHelper::InitialiseAnisotropicTestData(streamerSiteId, - streamerFOld); - - // Calculate what the value streamed to site streamedToSite should be. - lb::HydroVars streamerHydroVars(streamerFOld); - normalCollision->CalculatePreCollision(streamerHydroVars, streamedSite); - - normalCollision->Collide(&lbmParams, streamerHydroVars); - - // F_new should be equal to the value that was streamed - // from this other site in the same direction as we're - // streaming from. - REQUIRE(apprx(streamerHydroVars.GetFPostCollision()[streamedDirection]) - == streamedToFNew[streamedDirection]); - } else { - // The streamer index shows that no one has streamed to - // streamedToSite direction streamedDirection, therefore - // bounce back has happened in that site for that - // direction - - // Initialise streamedToSiteFOld with the original data - distribn_t streamerToSiteFOld[NUMVECTORS]; - LbTestsHelper::InitialiseAnisotropicTestData(streamedToSite, streamerToSiteFOld); - lb::HydroVars hydroVars(streamerToSiteFOld); - normalCollision->CalculatePreCollision(hydroVars, streamedSite); - - // Simulate post-collision using the collision operator. - normalCollision->Collide(&lbmParams, hydroVars); - - // After streaming FNew in a given direction must be f - // post-collision in the opposite direction following - // collision - INFO("Junk&Yang bounce-back equivalent: site " << streamedToSite - << " direction " << streamedDirection); - REQUIRE(apprx(streamedToFNew[streamedDirection]) - == hydroVars.GetFPostCollision()[oppDirection]); - } - } - } - } + SECTION("JunkYangEquivalentToBounceBack") { + // Initialise fOld in the lattice data. We choose values so + // that each site has an anisotropic distribution function, + // and that each site's function is distinguishable. + LbTestsHelper::InitialiseAnisotropicTestData(*latDat); + // Setting all the wall distances to 0.5 will make Junk&Yang + // behave like Simple Bounce Back + LbTestsHelper::SetWallAndIoletDistances(*latDat, 0.5); + + auto [firstBulkSite, pastLastBulkSite] = dom->GetMidDomainSiteRange(0); + auto wallSiteRange = dom->GetMidDomainSiteRange(1); + auto wallSitesCount = wallSiteRange.second - wallSiteRange.first; + // Check that the lattice has sites labeled as wall (otherwise + // this test is void) + REQUIRE(wallSitesCount == 24); + + // Mid-Fluid sites use simple collide and stream + BulkStreamer simpleCollideAndStream(initParams); + + simpleCollideAndStream.StreamAndCollide( + firstBulkSite, pastLastBulkSite, + &lbmParams, *latDat, *propertyCache + ); + + // Wall sites use Junk and Yang + initParams.siteRanges = {wallSiteRange}; + JunkYangFactory> junkYang(initParams); + + junkYang.StreamAndCollide( + wallSiteRange.first, wallSiteRange.second, + &lbmParams, *latDat, *propertyCache + ); + + junkYang.PostStep( + wallSiteRange.first, wallSiteRange.second, + &lbmParams, *latDat, *propertyCache + ); + + // Consider inlet/outlets and their walls as mid-fluid sites + simpleCollideAndStream.StreamAndCollide( + wallSiteRange.second, dom->GetLocalFluidSiteCount(), + &lbmParams, *latDat, *propertyCache + ); + + // Loop over the wall sites and check whether they got + // properly streamed on or bounced back depending on where + // they sit relative to the wall. We ignore mid-Fluid sites + // since StreamAndCollide was tested before. + for (site_t iSite = wallSiteRange.first; iSite < wallSiteRange.second; ++iSite) { + const auto streamedSite = latDat->GetSite(iSite); + distribn_t* streamedToFNew = latDat->GetFNew(NUMVECTORS * iSite); + + for (unsigned int iDir = 0; iDir < NUMVECTORS; ++iDir) { + unsigned oppDirection = LATTICE::INVERSEDIRECTIONS[iDir]; + + // Index of the site streaming to streamedToSite via + // direction streamedDirection + site_t streamerIndex = streamedSite.GetStreamedIndex (oppDirection); + + // Is streamerIndex a valid index? + if (streamerIndex >= 0 && + streamerIndex < (NUMVECTORS * dom->GetLocalFluidSiteCount())) { + // The streamer index is a valid index in the domain, + // therefore stream and collide has happened + site_t streamerSiteId = streamerIndex / NUMVECTORS; + + // Calculate streamerFOld at this site. + distribn_t streamerFOld[NUMVECTORS]; + LbTestsHelper::InitialiseAnisotropicTestData(streamerSiteId, + streamerFOld); + + // Calculate what the value streamed to site streamedToSite should be. + lb::HydroVars streamerHydroVars(streamerFOld); + normalCollision->CalculatePreCollision(streamerHydroVars, streamedSite); + + normalCollision->Collide(&lbmParams, streamerHydroVars); + + // F_new should be equal to the value that was streamed + // from this other site in the same direction as we're + // streaming from. + REQUIRE(apprx(streamerHydroVars.GetFPostCollision()[iDir]) + == streamedToFNew[iDir]); + } else { + // The streamer index shows that no one has streamed to + // streamedToSite direction streamedDirection, therefore + // bounce back has happened in that site for that + // direction + + // Initialise streamedToSiteFOld with the original data + distribn_t streamerToSiteFOld[NUMVECTORS]; + LbTestsHelper::InitialiseAnisotropicTestData(iSite, streamerToSiteFOld); + lb::HydroVars hydroVars(streamerToSiteFOld); + normalCollision->CalculatePreCollision(hydroVars, streamedSite); + + // Simulate post-collision using the collision operator. + normalCollision->Collide(&lbmParams, hydroVars); + + // After streaming FNew in a given direction must be f + // post-collision in the opposite direction following + // collision + INFO("Junk&Yang bounce-back equivalent: site " << iSite + << " direction " << iDir); + REQUIRE(apprx(streamedToFNew[iDir]) + == hydroVars.GetFPostCollision()[oppDirection]); + } + } + } + } SECTION("NashZerothOrderPressureIolet") { using N0P = StreamerTypeFactory, NashZerothOrderPressureLink>; diff --git a/Code/tests/lb/VirtualSiteIoletStreamerTests.cc b/Code/tests/lb/VirtualSiteIoletStreamerTests.cc index bfbaa689f..333885c7b 100644 --- a/Code/tests/lb/VirtualSiteIoletStreamerTests.cc +++ b/Code/tests/lb/VirtualSiteIoletStreamerTests.cc @@ -283,15 +283,7 @@ namespace hemelb::tests SECTION("TestStreamerInitialisation") { initParams.boundaryObject = &outletBoundary; // Set up the ranges to cover Mid 3 (pure outlet) and Mid 5 (outlet/wall) - initParams.siteRanges.resize(2); - initParams.siteRanges[0].first = dom->GetMidDomainCollisionCount(0) - + dom->GetMidDomainCollisionCount(1) + dom->GetMidDomainCollisionCount(2); - initParams.siteRanges[0].second = initParams.siteRanges[0].first - + dom->GetMidDomainCollisionCount(3); - initParams.siteRanges[1].first = initParams.siteRanges[0].second - + dom->GetMidDomainCollisionCount(4); - initParams.siteRanges[1].second = initParams.siteRanges[1].first - + dom->GetMidDomainCollisionCount(5); + initParams.siteRanges = {dom->GetMidDomainSiteRange(3), dom->GetMidDomainSiteRange(5)}; lb::VirtualSiteIolet outletStreamer(initParams); // All the sites at the outlet plane (x, y, 3) should be in the cache. @@ -331,32 +323,26 @@ namespace hemelb::tests InitialiseGradientHydroVars(); // Stream and collide - site_t offset = 0; - offset += dom->GetMidDomainCollisionCount(0); - offset += dom->GetMidDomainCollisionCount(1); - inletStreamer.StreamAndCollide(offset, - dom->GetMidDomainCollisionCount(2), - &lbmParams, - *latDat , - propertyCache); - offset += dom->GetMidDomainCollisionCount(2); - - outletStreamer.StreamAndCollide(offset, - dom->GetMidDomainCollisionCount(3), + inletStreamer.StreamAndCollide(dom->GetMidDomainSiteRange(2).first, + dom->GetMidDomainSiteRange(2).second, + &lbmParams, + *latDat , + propertyCache); + + outletStreamer.StreamAndCollide(dom->GetMidDomainSiteRange(3).first, + dom->GetMidDomainSiteRange(3).second, &lbmParams, *latDat, propertyCache); - offset += dom->GetMidDomainCollisionCount(3); - inletStreamer.StreamAndCollide(offset, - dom->GetMidDomainCollisionCount(4), + inletStreamer.StreamAndCollide(dom->GetMidDomainSiteRange(4).first, + dom->GetMidDomainSiteRange(4).second, &lbmParams, *latDat, propertyCache); - offset += dom->GetMidDomainCollisionCount(4); - outletStreamer.StreamAndCollide(offset, - dom->GetMidDomainCollisionCount(5), + outletStreamer.StreamAndCollide(dom->GetMidDomainSiteRange(5).first, + dom->GetMidDomainSiteRange(5).second, &lbmParams, *latDat, propertyCache); @@ -366,32 +352,26 @@ namespace hemelb::tests CheckAllHVUpdated(outletBoundary, 0); // Stream and collide - offset = 0; - offset += dom->GetMidDomainCollisionCount(0); - offset += dom->GetMidDomainCollisionCount(1); - inletStreamer.PostStep(offset, - dom->GetMidDomainCollisionCount(2), - &lbmParams, - latDat.get(), - propertyCache); - offset += dom->GetMidDomainCollisionCount(2); - - outletStreamer.PostStep(offset, - dom->GetMidDomainCollisionCount(3), + inletStreamer.PostStep(dom->GetMidDomainSiteRange(2).first, + dom->GetMidDomainSiteRange(2).second, + &lbmParams, + latDat.get(), + propertyCache); + + outletStreamer.PostStep(dom->GetMidDomainSiteRange(3).first, + dom->GetMidDomainSiteRange(3).second, &lbmParams, latDat.get(), propertyCache); - offset += dom->GetMidDomainCollisionCount(3); - inletStreamer.PostStep(offset, - dom->GetMidDomainCollisionCount(4), + inletStreamer.PostStep(dom->GetMidDomainSiteRange(4).first, + dom->GetMidDomainSiteRange(4).second, &lbmParams, latDat.get(), propertyCache); - offset += dom->GetMidDomainCollisionCount(4); - outletStreamer.PostStep(offset, - dom->GetMidDomainCollisionCount(5), + outletStreamer.PostStep(dom->GetMidDomainSiteRange(5).first, + dom->GetMidDomainSiteRange(5).second, &lbmParams, latDat.get(), propertyCache);