diff --git a/Code.v05-00/CMakeLists.txt b/Code.v05-00/CMakeLists.txt index 29d357d2e..ab286c11c 100644 --- a/Code.v05-00/CMakeLists.txt +++ b/Code.v05-00/CMakeLists.txt @@ -10,12 +10,13 @@ option(DEBUG "Enable debug mode" OFF) SET(RINGS 0) SET(OMP 1) option(BUILD_TEST ON) +option(ENABLE_PROFILING "Enable perf-friendly build (-g -O2)" OFF) if (CMAKE_BUILD_TYPE STREQUAL "Debug") include(CheckCXXCompilerFlag) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 -Wall -Wextra -Wconversion") #add_compile_options(-fno-omit-frame-pointer -fsanitize=bounds -fsanitize=address -fsanitize=undefined) - #add_link_options(-fsanitize=bounds -fsanitize=address -fsanitize=undefined) + #add_link_options(-fsanitize=bounds -fsanitize=address -fsanitize=undefined)s add_compile_options(-fno-omit-frame-pointer -fsanitize=bounds -fsanitize=undefined) add_link_options(-fsanitize=bounds -fsanitize=undefined) elseif (NOT CMAKE_BUILD_TYPE OR CMAKE_BUILD_TYPE STREQUAL "") @@ -34,6 +35,10 @@ elseif (NOT CMAKE_BUILD_TYPE OR CMAKE_BUILD_TYPE STREQUAL "") #-Michael #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ffast-math -fno-trapping-math -funsafe-math-optimizations") endif() +if(ENABLE_PROFILING) + message(STATUS "Profiling enabled: adding debug symbols and frame pointers") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O2 -fno-omit-frame-pointer") +endif() set(CMAKE_EXPORT_COMPILE_COMMANDS ON) execute_process(COMMAND git rev-parse --short HEAD @@ -131,4 +136,4 @@ add_subdirectory(tests) if (DEBUG) message(STATUS "DEBUG mode is enabled") -endif() +endif() \ No newline at end of file diff --git a/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp b/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp index 260350a0e..4060cd083 100644 --- a/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp +++ b/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp @@ -100,11 +100,31 @@ class LAGRIDPlumeModel { void updateDiffVecs(); void runTransport(double timestep); void remapAllVars(double remapTimestep, const std::vector>& mask, const VectorUtils::MaskInfo& maskInfo); - std::pair remapVariable(const VectorUtils::MaskInfo& maskInfo, const BufferInfo& buffers, const Vector_2D& phi, const std::vector>& mask); + std::pair remapVariable( + const VectorUtils::MaskInfo& maskInfo, + const BufferInfo& buffers, + const Vector_2D& phi, + const std::vector>& mask + ); double totalAirMass(); - Eigen::SparseMatrix createRegriddingWeightsSparse(const VectorUtils::MaskInfo& maskInfo, const BufferInfo& buffers, const std::vector>& mask, Vector_1D& xEdgesNew, Vector_1D& yEdgesNew, Vector_1D& xCoordsNew, Vector_1D& yCoordsNew); - Vector_2D applyWeights(const Eigen::SparseMatrix& weights, int nx_old, int ny_old, int nx_new, int ny_new, const Vector_2D& dataIn); + Eigen::SparseMatrix createRegriddingWeightsSparse( + const VectorUtils::MaskInfo& maskInfo, + const BufferInfo& buffers, + const std::vector>& mask, + Vector_1D& xEdgesNew, + Vector_1D& yEdgesNew, + Vector_1D& xCoordsNew, + Vector_1D& yCoordsNew + ); + Vector_2D applyWeights( + const Eigen::SparseMatrix& weights, + int nx_old, + int ny_old, + int nx_new, + int ny_new, + const Vector_2D& dataIn + ); void printVector2D(const std::string fieldName, const Vector_2D& dataIn); diff --git a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp index 6c6e06b20..050928cb5 100644 --- a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp +++ b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp @@ -5,8 +5,10 @@ #include "FVM_ANDS/BoundaryCondition.hpp" #include "FVM_ANDS_HelperFunctions.hpp" #include +#include namespace FVM_ANDS{ + using PointVariant = std::variant; // Separate the SOR solver for testing without having to build an AdvDiffSystem object void sor_solve(const Eigen::SparseMatrix &A, const Eigen::VectorXd &rhs, Eigen::VectorXd &phi, double omega = 1.0, double threshold = 1e-3, int n_iters = 3); @@ -41,7 +43,12 @@ namespace FVM_ANDS{ void sor_solve(double omega = 1.0, double threshold = 1e-3, int n_iters = 3){ FVM_ANDS::sor_solve(totalCoefMatrix_, rhs_, phi_, omega, threshold, n_iters); }; inline const Eigen::VectorXd& getRHS() const { return rhs_; } inline const Eigen::VectorXd& phi() const { return phi_; } - inline const std::vector>& points() const { return points_; } + inline const std::vector& points() const { return points_; } + inline const std::vector& bcCache() const { return bcCache_; } + inline const std::vector& directionCache() const { return directionCache_; } + inline const std::vector>& secondBoundaryCache() const { return secondBoundaryCache_; }; + inline const std::vector& corrCache() const { return corrCache_; } + inline const Eigen::SparseMatrix& getCoefMatrix() const { return totalCoefMatrix_; } inline void updatePhi(const Eigen::VectorXd& phi_new){ //Need to resize to account for grid changing in size. @@ -140,7 +147,11 @@ namespace FVM_ANDS{ Vector_1D bcVals_left_; Vector_1D bcVals_right_; Vector_1D bcVals_bot_; - std::vector> points_; + std::vector points_; + std::vector bcCache_; // Cache bcType() + std::vector directionCache_; // Cache bcDirection() + std::vector> secondBoundaryCache_; // Cache secondBoundaryConds() + std::vector corrCache_; // Cache corrPoint() Eigen::SparseMatrix totalCoefMatrix_; Eigen::VectorXd rhs_; Eigen::VectorXd phi_; @@ -151,6 +162,12 @@ namespace FVM_ANDS{ void buildPointList(); void buildAdvectionCoeffs(int i, double& coeff_C, double& coeff_N, double& coeff_S, double& coeff_E, double& coeff_W); void updateGhostNodes(); + + // Helper to create points + template + void addPoint(int idx, Args&&... args) { + points_[idx] = T{std::forward(args)...}; + } inline bool isValidPointID(int idx) const { return (idx >= 0 && idx < phi_.rows()); @@ -306,14 +323,14 @@ namespace FVM_ANDS{ return std::max(0.0, std::min(r, 1.0)); } inline int neighbor_point(FaceDirection direction, int pointID) const noexcept{ - Point* point = points_[pointID].get(); - if(point->bcType() == BoundaryConditionFlag::INTERIOR) return neighbor_point_interior(direction, pointID); + Point point = points_[pointID]; + if(point.bcType() == BoundaryConditionFlag::INTERIOR) return neighbor_point_interior(direction, pointID); - if(point->bcDirection() == direction){ - return point->corrPoint(); + if(point.bcDirection() == direction){ + return point.corrPoint(); } - else if (point->secondBoundaryConds() && point->secondBoundaryConds().value().direction == direction){ - return point->secondBoundaryConds().value().corrPoint; + else if (point.secondBoundaryConds() && point.secondBoundaryConds().value().direction == direction){ + return point.secondBoundaryConds().value().corrPoint; } return neighbor_point_interior(direction, pointID); } diff --git a/Code.v05-00/include/FVM_ANDS/BoundaryCondition.hpp b/Code.v05-00/include/FVM_ANDS/BoundaryCondition.hpp index fc6d0b85e..8027c0557 100644 --- a/Code.v05-00/include/FVM_ANDS/BoundaryCondition.hpp +++ b/Code.v05-00/include/FVM_ANDS/BoundaryCondition.hpp @@ -44,7 +44,7 @@ namespace FVM_ANDS{ inline virtual void setBCType(BoundaryConditionFlag bcType){ bc_.bcType = bcType; } - constexpr virtual bool isGhost(){ + constexpr virtual bool isGhost() const { return false; } inline virtual void setSecondaryBC(BoundaryCondDescription bc) { @@ -60,7 +60,7 @@ namespace FVM_ANDS{ public: GhostPoint() = delete; GhostPoint(BoundaryCondDescription bc); - constexpr bool isGhost() override { + constexpr bool isGhost() const override { return true; } inline void setBCType(BoundaryConditionFlag bcType) override { diff --git a/Code.v05-00/include/FVM_ANDS/FVM_Solver.hpp b/Code.v05-00/include/FVM_ANDS/FVM_Solver.hpp index 27c78ca65..9da529274 100644 --- a/Code.v05-00/include/FVM_ANDS/FVM_Solver.hpp +++ b/Code.v05-00/include/FVM_ANDS/FVM_Solver.hpp @@ -53,7 +53,7 @@ namespace FVM_ANDS{ inline const Eigen::SparseMatrix& coefMatrix(){ return advDiffSys_.getCoefMatrix(); } - inline const std::vector>& points(){ + inline const std::vector& points(){ return advDiffSys_.points(); } const Eigen::VectorXd& calcRHS(){ diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index 90edda0d2..861fede97 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -1,6 +1,7 @@ #include #include #include +#include namespace FVM_ANDS{ AdvDiffSystem::AdvDiffSystem(const AdvDiffParams& params, const Vector_1D xCoords, const Vector_1D yCoords, const BoundaryConditions& bc, const Eigen::VectorXd& phi_init, vecFormat format) : @@ -36,12 +37,24 @@ namespace FVM_ANDS{ Dv_vec_.resize(nInteriorPoints_); rhs_.resize(nTotalPoints_); phi_.resize(nTotalPoints_); - points_.reserve(nTotalPoints_); + points_.resize(nTotalPoints_); + bcCache_.resize(nTotalPoints_); + directionCache_.resize(nTotalPoints_); + secondBoundaryCache_.resize(nTotalPoints_); + corrCache_.resize(nTotalPoints_); deferredCorr_.resize(nInteriorPoints_); deferredCorr_.setZero(); source_.resize(nInteriorPoints_); source_.setZero(); - std::generate_n(std::back_inserter(points_), nTotalPoints_, [] { return std::make_unique(); }); + std::generate_n(std::back_inserter(points_), nTotalPoints_, [] { + return Point(); // Construct directly in vector + }); + for (size_t i = 0; i < points_.size(); ++i) { + bcCache_[i] = points_[i].bcType(); + directionCache_[i] = points_[i].bcDirection(); + secondBoundaryCache_[i] = points_[i].secondBoundaryConds(); + corrCache_[i] = points_[i].corrPoint(); + } totalCoefMatrix_.resize(nTotalPoints_, nTotalPoints_); updateDiffusion(params.Dh, params.Dv); @@ -69,7 +82,7 @@ namespace FVM_ANDS{ if(j == ny_ - 1){ BoundaryCondDescription bc_ghost = BoundaryCondDescription(bcType_top_, FaceDirection::NORTH, bcVals_top_[i], vector_idx); int corrGhostPoint = nInteriorPoints_ + i; - points_[corrGhostPoint] = std::make_unique(bc_ghost); + points_[corrGhostPoint] = GhostPoint{bc_ghost}; BoundaryCondDescription bc_top = BoundaryCondDescription(bcType_top_, FaceDirection::NORTH, bcVals_top_[i], corrGhostPoint); @@ -79,27 +92,27 @@ namespace FVM_ANDS{ int corrGhostPoint2 = nInteriorPoints_ + nx_ + j; BoundaryCondDescription bc_left = BoundaryCondDescription(bcType_left_, FaceDirection::WEST, bcVals_left_[j], corrGhostPoint2); BoundaryCondDescription bc_ghost2 = BoundaryCondDescription(bcType_left_, FaceDirection::WEST, bcVals_left_[j], vector_idx); - points_[vector_idx] = std::make_unique(bc_top, bc_left); - points_[corrGhostPoint2] = std::make_unique(bc_ghost2); + points_[vector_idx] = IntBoundPoint(bc_top, bc_left); + points_[corrGhostPoint2] = GhostPoint(bc_ghost2); } //top right else if(i == nx_ - 1){ int corrGhostPoint2 = nInteriorPoints_ + nx_ + ny_ + j; BoundaryCondDescription bc_right = BoundaryCondDescription(bcType_right_, FaceDirection::EAST, bcVals_right_[j], corrGhostPoint2); BoundaryCondDescription bc_ghost2 = BoundaryCondDescription(bcType_right_, FaceDirection::EAST, bcVals_right_[j], vector_idx); - points_[vector_idx] = std::make_unique(bc_top, bc_right); - points_[corrGhostPoint2] = std::make_unique(bc_ghost2); + points_[vector_idx] = IntBoundPoint(bc_top, bc_right); + points_[corrGhostPoint2] = GhostPoint(bc_ghost2); } //if not corner boundary make normal interior node. else{ - points_[vector_idx] = std::make_unique(bc_top); + points_[vector_idx] = IntBoundPoint(bc_top); } } //bottom else if (j == 0){ int corrGhostPoint = nInteriorPoints_ + nx_ + 2*ny_ + i; BoundaryCondDescription bc_ghost = BoundaryCondDescription(bcType_bot_, FaceDirection::SOUTH, bcVals_bot_[i], vector_idx); - points_[corrGhostPoint] = std::make_unique(bc_ghost); + points_[corrGhostPoint] = GhostPoint(bc_ghost); BoundaryCondDescription bc_bot = BoundaryCondDescription(bcType_bot_, FaceDirection::SOUTH, bcVals_bot_[i], corrGhostPoint); @@ -111,8 +124,8 @@ namespace FVM_ANDS{ BoundaryCondDescription bc_left = BoundaryCondDescription(bcType_left_, FaceDirection::WEST, bcVals_left_[j], corrGhostPoint2); BoundaryCondDescription bc_ghost2 = BoundaryCondDescription(bcType_left_, FaceDirection::WEST, bcVals_left_[j], vector_idx); - points_[vector_idx] = std::make_unique(bc_bot, bc_left); - points_[corrGhostPoint2] = std::make_unique(bc_ghost2); + points_[vector_idx] = IntBoundPoint(bc_bot, bc_left); + points_[corrGhostPoint2] = GhostPoint(bc_ghost2); } //bot right else if(i == nx_ - 1){ @@ -120,13 +133,13 @@ namespace FVM_ANDS{ BoundaryCondDescription bc_right = BoundaryCondDescription(bcType_right_, FaceDirection::EAST, bcVals_right_[j], corrGhostPoint2); BoundaryCondDescription bc_ghost2 = BoundaryCondDescription(bcType_right_, FaceDirection::EAST, bcVals_right_[j], vector_idx); - points_[vector_idx] = std::make_unique(bc_bot, bc_right); - points_[corrGhostPoint2] = std::make_unique(bc_ghost2); + points_[vector_idx] = IntBoundPoint(bc_bot, bc_right); + points_[corrGhostPoint2] = GhostPoint(bc_ghost2); } //if not corner boundary make normal interior node else { - points_[vector_idx] = std::make_unique(bc_bot); + points_[vector_idx] = IntBoundPoint(bc_bot); } } @@ -136,8 +149,8 @@ namespace FVM_ANDS{ BoundaryCondDescription bc_left = BoundaryCondDescription(bcType_left_, FaceDirection::WEST, bcVals_left_[j], corrGhostPoint); BoundaryCondDescription bc_ghost = BoundaryCondDescription(bcType_left_, FaceDirection::WEST, bcVals_left_[j], vector_idx); - points_[vector_idx] = std::make_unique(bc_left); - points_[corrGhostPoint] = std::make_unique(bc_ghost); + points_[vector_idx] = IntBoundPoint(bc_left); + points_[corrGhostPoint] = GhostPoint(bc_ghost); } //right else if (i == nx_ - 1){ @@ -145,11 +158,11 @@ namespace FVM_ANDS{ BoundaryCondDescription bc_right = BoundaryCondDescription(bcType_right_, FaceDirection::EAST, bcVals_right_[j], corrGhostPoint); BoundaryCondDescription bc_ghost = BoundaryCondDescription(bcType_right_, FaceDirection::EAST, bcVals_right_[j], vector_idx); - points_[vector_idx] = std::make_unique(bc_right); - points_[corrGhostPoint] = std::make_unique(bc_ghost); + points_[vector_idx] = IntBoundPoint(bc_right); + points_[corrGhostPoint] = GhostPoint(bc_ghost); } else { - points_[vector_idx] = std::make_unique(BoundaryConditionFlag::INTERIOR); + points_[vector_idx] = Point(BoundaryConditionFlag::INTERIOR); } } } @@ -163,13 +176,13 @@ namespace FVM_ANDS{ auto start = std::chrono::high_resolution_clock::now(); for(int i = 0; i < nTotalPoints_; i++){ - if(points_[i]->isGhost()){ - switch(points_[i]->bcType()){ + if(points_[i].isGhost()){ + switch(points_[i].bcType()){ case BoundaryConditionFlag::DIRICHLET_GHOSTPOINT:{ // (phi_int + phi_ghost) / 2 = phi_boundary // if inhomog, the bc value will appear in the rhs. tripletList.emplace_back(i, i, 0.5); - tripletList.emplace_back(i, points_[i]->corrPoint(), 0.5); + tripletList.emplace_back(i, points_[i].corrPoint(), 0.5); break; } case BoundaryConditionFlag::PERIODIC_GHOSTPOINT:{ @@ -220,16 +233,16 @@ namespace FVM_ANDS{ //Therefore, that term goes to the RHS and the contribution of that face to the coeffs goes to 0. bool isNorthBoundary = 0, isWestBoundary = 0, isEastBoundary = 0, isSouthBoundary = 0; - if(points_[i]->bcType() != BoundaryConditionFlag::INTERIOR){ - isNorthBoundary = points_[i]->bcDirection() == FaceDirection::NORTH; - isSouthBoundary = points_[i]->bcDirection() == FaceDirection::SOUTH; + if(points_[i].bcType() != BoundaryConditionFlag::INTERIOR){ + isNorthBoundary = points_[i].bcDirection() == FaceDirection::NORTH; + isSouthBoundary = points_[i].bcDirection() == FaceDirection::SOUTH; //Corner cases... - bool secondaryWestBound = (points_[i]->secondBoundaryConds() && points_[i]->secondBoundaryConds().value().direction == FaceDirection::WEST); - bool secondaryEastBound = (points_[i]->secondBoundaryConds() && points_[i]->secondBoundaryConds().value().direction == FaceDirection::EAST); + bool secondaryWestBound = (points_[i].secondBoundaryConds() && points_[i].secondBoundaryConds().value().direction == FaceDirection::WEST); + bool secondaryEastBound = (points_[i].secondBoundaryConds() && points_[i].secondBoundaryConds().value().direction == FaceDirection::EAST); - isWestBoundary = (points_[i]->bcDirection() == FaceDirection::WEST || secondaryWestBound); - isEastBoundary = (points_[i]->bcDirection() == FaceDirection::EAST || secondaryEastBound); + isWestBoundary = (points_[i].bcDirection() == FaceDirection::WEST || secondaryWestBound); + isEastBoundary = (points_[i].bcDirection() == FaceDirection::EAST || secondaryEastBound); } int idx_E = neighbor_point(FaceDirection::EAST, i); @@ -277,11 +290,11 @@ namespace FVM_ANDS{ const Eigen::VectorXd& AdvDiffSystem::calcRHS(){ for(int i = 0; i < nTotalPoints_; i++){ - if(points_[i]->isGhost()){ - switch(points_[i]->bcType()){ + if(points_[i].isGhost()){ + switch(points_[i].bcType()){ case BoundaryConditionFlag::DIRICHLET_GHOSTPOINT:{ // Equation: (phi_int + phi_ghost) / 2 = phi_boundary - rhs_[i] = points_[i]->bcVal(); + rhs_[i] = points_[i].bcVal(); break; } case BoundaryConditionFlag::PERIODIC_GHOSTPOINT:{ @@ -297,31 +310,31 @@ namespace FVM_ANDS{ continue; } - switch(points_[i]->bcType()){ + switch(points_[i].bcType()){ case BoundaryConditionFlag::INTERIOR:{ rhs_[i] = phi_[i] + deferredCorr_[i] + source_[i]*dt_; break; } case BoundaryConditionFlag::DIRICHLET_INT_BPOINT:{ rhs_[i] = phi_[i] + deferredCorr_[i] + source_[i]*dt_; - switch(points_[i]->bcDirection()){ + switch(points_[i].bcDirection()){ case FaceDirection::NORTH: - rhs_[i] -= v_vec_[i] * dt_ / dy_ * points_[i]->bcVal(); + rhs_[i] -= v_vec_[i] * dt_ / dy_ * points_[i].bcVal(); break; case FaceDirection::SOUTH: - rhs_[i] += v_vec_[i] * dt_ / dy_ * points_[i]->bcVal(); + rhs_[i] += v_vec_[i] * dt_ / dy_ * points_[i].bcVal(); break; case FaceDirection::EAST: - rhs_[i] -= u_vec_[i] * dt_ / dx_ * points_[i]->bcVal(); + rhs_[i] -= u_vec_[i] * dt_ / dx_ * points_[i].bcVal(); break; case FaceDirection::WEST: - rhs_[i] += u_vec_[i] * dt_ / dx_ * points_[i]->bcVal(); + rhs_[i] += u_vec_[i] * dt_ / dx_ * points_[i].bcVal(); break; case FaceDirection::ERROR: throw std::runtime_error("Invalid FaceDirection in Dirichlet boundary condition"); } - if (!points_[i]->secondBoundaryConds()) break; - BoundaryCondDescription bc_2 = points_[i]->secondBoundaryConds().value(); + if (!points_[i].secondBoundaryConds()) break; + BoundaryCondDescription bc_2 = points_[i].secondBoundaryConds().value(); switch(bc_2.direction){ case FaceDirection::EAST: rhs_[i] -= u_vec_[i] * dt_ / dx_ * bc_2.bcVal; @@ -354,8 +367,8 @@ namespace FVM_ANDS{ int bPointID_top = twoDIdx_to_vecIdx(i, ny_ - 1, nx_, ny_, format_); switch(bcType_top_){ case BoundaryConditionFlag::DIRICHLET_INT_BPOINT: { - int ghostPointID = points_[bPointID_top]->corrPoint(); - phi_[ghostPointID] = 2 * points_[bPointID_top]->bcVal() - phi_[bPointID_top]; + int ghostPointID = points_[bPointID_top].corrPoint(); + phi_[ghostPointID] = 2 * points_[bPointID_top].bcVal() - phi_[bPointID_top]; break; } default: { @@ -367,8 +380,8 @@ namespace FVM_ANDS{ int bPointID_bot = twoDIdx_to_vecIdx(i, 0, nx_, ny_, format_); switch(bcType_bot_){ case BoundaryConditionFlag::DIRICHLET_INT_BPOINT: { - int ghostPointID = points_[bPointID_bot]->corrPoint(); - phi_[ghostPointID] = 2 * points_[bPointID_bot]->bcVal() - phi_[bPointID_bot]; + int ghostPointID = points_[bPointID_bot].corrPoint(); + phi_[ghostPointID] = 2 * points_[bPointID_bot].bcVal() - phi_[bPointID_bot]; break; } default: { @@ -382,8 +395,8 @@ namespace FVM_ANDS{ //corner cases if(j == 0 || j == ny_ - 1){ int bPointID_cornerLeft = twoDIdx_to_vecIdx(0, j, nx_, ny_, format_); - int ghostPointID = points_[bPointID_cornerLeft]->secondBoundaryConds().value().corrPoint; - double bcVal = points_[bPointID_cornerLeft]->secondBoundaryConds().value().bcVal; + int ghostPointID = points_[bPointID_cornerLeft].secondBoundaryConds().value().corrPoint; + double bcVal = points_[bPointID_cornerLeft].secondBoundaryConds().value().bcVal; switch(bcType_left_){ case BoundaryConditionFlag::DIRICHLET_INT_BPOINT: { phi_[ghostPointID] = 2 * bcVal - phi_[bPointID_cornerLeft]; @@ -395,8 +408,8 @@ namespace FVM_ANDS{ } int bPointID_cornerRight = twoDIdx_to_vecIdx(nx_ - 1, j, nx_, ny_, format_); - ghostPointID = points_[bPointID_cornerRight]->secondBoundaryConds().value().corrPoint; - bcVal = points_[bPointID_cornerRight]->secondBoundaryConds().value().bcVal; + ghostPointID = points_[bPointID_cornerRight].secondBoundaryConds().value().corrPoint; + bcVal = points_[bPointID_cornerRight].secondBoundaryConds().value().bcVal; switch(bcType_right_){ case BoundaryConditionFlag::DIRICHLET_INT_BPOINT: { @@ -413,8 +426,8 @@ namespace FVM_ANDS{ int bPointID_left = twoDIdx_to_vecIdx(0, j, nx_, ny_, format_); switch(bcType_left_){ case BoundaryConditionFlag::DIRICHLET_INT_BPOINT: { - int ghostPointID = points_[bPointID_left]->corrPoint(); - phi_[ghostPointID] = 2 * points_[bPointID_left]->bcVal() - phi_[bPointID_left]; + int ghostPointID = points_[bPointID_left].corrPoint(); + phi_[ghostPointID] = 2 * points_[bPointID_left].bcVal() - phi_[bPointID_left]; break; } default: { @@ -425,8 +438,8 @@ namespace FVM_ANDS{ int bPointID_right = twoDIdx_to_vecIdx(nx_ - 1, j, nx_, ny_, format_); switch(bcType_right_){ case BoundaryConditionFlag::DIRICHLET_INT_BPOINT: { - int ghostPointID = points_[bPointID_right]->corrPoint(); - phi_[ghostPointID] = 2 * points_[bPointID_right]->bcVal() - phi_[bPointID_right]; + int ghostPointID = points_[bPointID_right].corrPoint(); + phi_[ghostPointID] = 2 * points_[bPointID_right].bcVal() - phi_[bPointID_right]; break; } default: { @@ -450,54 +463,54 @@ namespace FVM_ANDS{ int currIdx = nInteriorPoints_; //top for(int i = 0; i < nx_; i++){ - int corrPointID = points_[currIdx]->corrPoint(); - points_[currIdx]->setBCType(bcType_top_); - points_[currIdx]->setBCVal(bcVals_top_[i]); - points_[corrPointID]->setBCType(bcType_top_); - points_[corrPointID]->setBCVal(bcVals_top_[i]); + int corrPointID = points_[currIdx].corrPoint(); + points_[currIdx].setBCType(bcType_top_); + points_[currIdx].setBCVal(bcVals_top_[i]); + points_[corrPointID].setBCType(bcType_top_); + points_[corrPointID].setBCVal(bcVals_top_[i]); currIdx++; } //left for(int i = 0; i < ny_; i++){ - int corrPointID = points_[currIdx]->corrPoint(); + int corrPointID = points_[currIdx].corrPoint(); if(i == 0 || i == ny_ - 1){ - points_[currIdx]->setBCType(bcType_left_); - points_[currIdx]->setBCVal(bcVals_left_[i]); + points_[currIdx].setBCType(bcType_left_); + points_[currIdx].setBCVal(bcVals_left_[i]); BoundaryCondDescription bc(bcType_left_, FaceDirection::WEST, bcVals_left_[i], currIdx); - points_[corrPointID]->setSecondaryBC(bc); + points_[corrPointID].setSecondaryBC(bc); currIdx++; continue; } - points_[currIdx]->setBCType(bcType_left_); - points_[currIdx]->setBCVal(bcVals_left_[i]); - points_[corrPointID]->setBCType(bcType_left_); - points_[corrPointID]->setBCVal(bcVals_left_[i]); + points_[currIdx].setBCType(bcType_left_); + points_[currIdx].setBCVal(bcVals_left_[i]); + points_[corrPointID].setBCType(bcType_left_); + points_[corrPointID].setBCVal(bcVals_left_[i]); currIdx++; } //right for(int i = 0; i < ny_; i++){ - int corrPointID = points_[currIdx]->corrPoint(); + int corrPointID = points_[currIdx].corrPoint(); if(i == 0 || i == ny_ - 1){ - points_[currIdx]->setBCType(bcType_right_); - points_[currIdx]->setBCVal(bcVals_right_[i]); + points_[currIdx].setBCType(bcType_right_); + points_[currIdx].setBCVal(bcVals_right_[i]); BoundaryCondDescription bc(bcType_right_, FaceDirection::EAST, bcVals_right_[i], currIdx); - points_[corrPointID]->setSecondaryBC(bc); + points_[corrPointID].setSecondaryBC(bc); currIdx++; continue; } - points_[currIdx]->setBCType(bcType_right_); - points_[currIdx]->setBCVal(bcVals_right_[i]); - points_[corrPointID]->setBCType(bcType_right_); - points_[corrPointID]->setBCVal(bcVals_right_[i]); + points_[currIdx].setBCType(bcType_right_); + points_[currIdx].setBCVal(bcVals_right_[i]); + points_[corrPointID].setBCType(bcType_right_); + points_[corrPointID].setBCVal(bcVals_right_[i]); currIdx++; } //bot for(int i = 0; i < nx_; i++){ - int corrPointID = points_[currIdx]->corrPoint(); - points_[currIdx]->setBCType(bcType_bot_); - points_[currIdx]->setBCVal(bcVals_bot_[i]); - points_[corrPointID]->setBCType(bcType_bot_); - points_[corrPointID]->setBCVal(bcVals_bot_[i]); + int corrPointID = points_[currIdx].corrPoint(); + points_[currIdx].setBCType(bcType_bot_); + points_[currIdx].setBCVal(bcVals_bot_[i]); + points_[corrPointID].setBCType(bcType_bot_); + points_[corrPointID].setBCVal(bcVals_bot_[i]); currIdx++; } applyBoundaryCondition(); //need this to calculate minmod function at some timestep. @@ -522,24 +535,44 @@ namespace FVM_ANDS{ //commenting out this results in ~30% speedup //The calls involving the optional are maybe 1/3 of the cost. Maybe something to look at later. - if(points_[i]->bcType() != BoundaryConditionFlag::INTERIOR){ - Point* point = points_[i].get(); - FaceDirection direction = point->bcDirection(); + if(bcCache_[i] != BoundaryConditionFlag::INTERIOR){ + // const Point& point = points_[i]; + // FaceDirection direction = point->bcDirection(); + // isNorthBoundary = direction == FaceDirection::NORTH; + // isSouthBoundary = direction == FaceDirection::SOUTH; + + FaceDirection direction = directionCache_[i]; isNorthBoundary = direction == FaceDirection::NORTH; isSouthBoundary = direction == FaceDirection::SOUTH; + bool secondaryWestBound = false; + bool secondaryEastBound = false; + + std::optional secondBC = secondBoundaryCache_[i]; + //Corner cases... - bool secondaryWestBound = (point->secondBoundaryConds() && point->secondBoundaryConds()->direction == FaceDirection::WEST); - bool secondaryEastBound = (point->secondBoundaryConds() && point->secondBoundaryConds()->direction == FaceDirection::EAST); + if (secondBC) { + secondaryWestBound = (secondBC->direction == FaceDirection::WEST); + secondaryEastBound = (secondBC->direction == FaceDirection::EAST); + } isWestBoundary = (direction == FaceDirection::WEST || secondaryWestBound); isEastBoundary = (direction == FaceDirection::EAST || secondaryEastBound); //only call this lookup function on boundary nodes which are inconsequential in number - idx_N = isNorthBoundary? point->corrPoint() : idx_N; - idx_S = isSouthBoundary? point->corrPoint() : idx_S; - idx_E = isEastBoundary? (secondaryEastBound ? point->secondBoundaryConds()->corrPoint : point->corrPoint()) : idx_E; - idx_W = isWestBoundary? (secondaryEastBound ? point->secondBoundaryConds()->corrPoint : point->corrPoint()) : idx_W; + idx_N = isNorthBoundary? corrCache_[i] : idx_N; + idx_S = isSouthBoundary? corrCache_[i] : idx_S; + if (isEastBoundary && secondaryEastBound && secondBC) { + idx_E = secondBC->corrPoint; + } else if (isEastBoundary) { + idx_E = corrCache_[i]; + } + + if (isWestBoundary && secondaryWestBound && secondBC) { + idx_W = secondBC->corrPoint; + } else if (isWestBoundary) { + idx_W = corrCache_[i]; + } } //When you declare these vars (inside or outside loop) has 0 impact) //takes ~ 6 out of 18 ns on background var calcs @@ -570,7 +603,7 @@ namespace FVM_ANDS{ //Using only first order upwind can result in a ~40% speedup of the total advection calc. //So... there is significantly more cost from actually doing the calculation than from branching. if(isNorthBoundary){ - phi_N = points_[i]->bcVal(); + phi_N = points_[i].bcVal(); } else if (v_local >= 0){ phi_N = phi_[i] + 0.5 * minmod_N_vPos(i) * (phi_[idx_N] - phi_[i]); @@ -579,7 +612,7 @@ namespace FVM_ANDS{ phi_N = phi_[idx_N] + 0.5 * minmod_N_vNeg(i) * (phi_[i] - phi_[idx_N]); } if(isSouthBoundary){ - phi_S = points_[i]->bcVal(); + phi_S = points_[i].bcVal(); } else if (v_local >= 0){ phi_S = phi_[idx_S] + 0.5 * minmod_S_vPos(i) * (phi_[i] - phi_[idx_S]); @@ -589,7 +622,7 @@ namespace FVM_ANDS{ } if(isWestBoundary){ - phi_W = secondaryWestBound ? points_[i]->secondBoundaryConds().value().bcVal : points_[i]->bcVal(); + phi_W = secondaryWestBound ? points_[i].secondBoundaryConds().value().bcVal : points_[i].bcVal(); } else if (u_local >= 0){ phi_W = phi_[idx_W] + 0.5 * minmod_W_vPos(i) * (phi_[i] - phi_[idx_W]); @@ -599,7 +632,7 @@ namespace FVM_ANDS{ } if(isEastBoundary){ - phi_E = secondaryEastBound ? points_[i]->secondBoundaryConds().value().bcVal : points_[i]->bcVal(); + phi_E = secondaryEastBound ? points_[i].secondBoundaryConds().value().bcVal : points_[i].bcVal(); } else if (u_local >= 0){ phi_E = phi_[i] + 0.5 * minmod_E_vPos(i) * (phi_[idx_E] - phi_[i]); diff --git a/Code.v05-00/tests/test_adv_diff_solver.cpp b/Code.v05-00/tests/test_adv_diff_solver.cpp index 286fdee55..219b6cbcd 100644 --- a/Code.v05-00/tests/test_adv_diff_solver.cpp +++ b/Code.v05-00/tests/test_adv_diff_solver.cpp @@ -159,7 +159,7 @@ namespace FVM_ANDS{ solver.buildCoeffMatrix(); auto mat = solver.coefMatrix(); - const std::vector>& points = solver.points(); + const std::vector& points = solver.points(); SECTION("Basic Sanity Checks"){ REQUIRE(mat.cols() == 26); @@ -168,8 +168,8 @@ namespace FVM_ANDS{ int numGhost = 0; int numBound = 0; for(int i = 0; i < 26; i++){ - if(points[i]->isGhost()) {numGhost++;} - else if (points[i]->bcType() != BoundaryConditionFlag::INTERIOR) numBound++; + if(points[i].isGhost()) {numGhost++;} + else if (points[i].bcType() != BoundaryConditionFlag::INTERIOR) numBound++; } REQUIRE(numGhost == 14); REQUIRE(numBound == 10); diff --git a/examples/issl_rhi140/input.yaml b/examples/issl_rhi140/input.yaml index b259c975e..8eafca27f 100644 --- a/examples/issl_rhi140/input.yaml +++ b/examples/issl_rhi140/input.yaml @@ -2,7 +2,7 @@ SIMULATION MENU: # Only one of parameter sweep or MC simulation can be on. # Parameter sweep lets you specify an arbitrary number of custom values for each parameter; Monte Carlo simulation is self-explanatory. # At the moment, you cannot mix and match MC sim and param sweep on each individual parameter. - OpenMP Num Threads (positive int): 8 + OpenMP Num Threads (positive int): 1 PARAM SWEEP SUBMENU: Parameter sweep (T/F): T #-OR--------------- @@ -32,7 +32,7 @@ SIMULATION MENU: Run box model (T/F): F netCDF filename format (string): APCEMM_BOX_CASE_* RANDOM NUMBER GENERATION SUBMENU: - Force seed value (T/F): F + Force seed value (T/F): T Seed value (positive int): 0 EPM type (original/external/new): original