From 55567a7400cadb2db97ca881765bbbab50c86238 Mon Sep 17 00:00:00 2001 From: Coco Yeung Date: Mon, 8 Dec 2025 16:28:36 +0000 Subject: [PATCH 1/8] Add profiling flags --- Code.v05-00/CMakeLists.txt | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Code.v05-00/CMakeLists.txt b/Code.v05-00/CMakeLists.txt index 29d357d2..ab286c11 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 From 04d337b16ccf99d3015f3c276e642554b1b3494c Mon Sep 17 00:00:00 2001 From: Coco Yeung Date: Mon, 8 Dec 2025 17:12:01 +0000 Subject: [PATCH 2/8] standardise seed and set threads to 1 --- examples/issl_rhi140/input.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/issl_rhi140/input.yaml b/examples/issl_rhi140/input.yaml index b259c975..8eafca27 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 From fce950fed596184f976073349607861801ae4f24 Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Tue, 30 Dec 2025 16:30:38 +0800 Subject: [PATCH 3/8] Getting rid of unique_ptr --- Code.v05-00/include/Core/LAGRIDPlumeModel.hpp | 26 +++++- .../include/FVM_ANDS/AdvDiffSystem.hpp | 21 ++++- Code.v05-00/src/Core/LAGRIDPlumeModel.cpp | 2 +- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 80 ++++++++++++------- 4 files changed, 94 insertions(+), 35 deletions(-) diff --git a/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp b/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp index 260350a0..b4968c42 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, V + ector_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 6c6e06b2..2b3d4cc0 100644 --- a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp +++ b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp @@ -5,6 +5,8 @@ #include "FVM_ANDS/BoundaryCondition.hpp" #include "FVM_ANDS_HelperFunctions.hpp" #include +#include +using PointVariant = std::variant; namespace FVM_ANDS{ @@ -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> secondBounaryCache_; // 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()); diff --git a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp index bcb0db4e..eb36eca1 100644 --- a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp +++ b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp @@ -367,7 +367,7 @@ void LAGRIDPlumeModel::updateDiffVecs() { } } } -void LAGRIDPlumeModel::runTransport(double timestep) { +void LAGRIDPlumeModel::runTrans,port(double timestep) { //Update the zero bc to reflect grid size changes auto ZERO_BC = FVM_ANDS::bcFrom2DVector(iceAerosol_.getPDF()[0], true); diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index 90edda0d..1c3a51bf 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -36,12 +36,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) { + bcTypeCache_[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 +81,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 +91,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 +123,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 +132,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 +148,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 +157,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); } } } @@ -505,6 +517,12 @@ namespace FVM_ANDS{ Eigen::VectorXd AdvDiffSystem::forwardEulerAdvection(bool operatorSplit, bool parallelAdvection) const noexcept{ Eigen::VectorXd soln(nTotalPoints_); + // cache results outside of hotloop + for (const auto& point : points_) { + bcTypes_.push_back(std::visit([](auto&& p) { return p.bcType(); }, point)); + directions_.push_back(std::visit([](auto&& p) { return p.bcDirection(); }, point)); + } + // double avgBackgroundCalcTime = 0; //Explicit Time-Stepping #pragma omp parallel for \ @@ -522,24 +540,28 @@ 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 = directions_[i]; isNorthBoundary = direction == FaceDirection::NORTH; isSouthBoundary = direction == FaceDirection::SOUTH; - + //Corner cases... - bool secondaryWestBound = (point->secondBoundaryConds() && point->secondBoundaryConds()->direction == FaceDirection::WEST); - bool secondaryEastBound = (point->secondBoundaryConds() && point->secondBoundaryConds()->direction == FaceDirection::EAST); + bool secondaryWestBound = (secondBoundaryCache_[i]->direction == FaceDirection::WEST); + bool secondaryEastBound = (secondBoundaryCache_[i]->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; + idx_E = isEastBoundary? (secondaryEastBound ? secondBoundaryCache_[i]->corrPoint : corrCache_[i]) : idx_E; + idx_W = isWestBoundary? (secondaryEastBound ? secondBoundaryCache_[i]->corrPoint : corrCache_[i]) : idx_W; } //When you declare these vars (inside or outside loop) has 0 impact) //takes ~ 6 out of 18 ns on background var calcs From 71290dd4f02d7ecbaab656f8e736b67e3c401e38 Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Tue, 30 Dec 2025 16:41:21 +0800 Subject: [PATCH 4/8] remove -> --- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 142 ++++++++++----------- 1 file changed, 71 insertions(+), 71 deletions(-) diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index 1c3a51bf..058012a6 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -49,10 +49,10 @@ namespace FVM_ANDS{ return Point(); // Construct directly in vector }); for (size_t i = 0; i < points_.size(); ++i) { - bcTypeCache_[i] = points_[i]->bcType(); - directionCache_[i] = points_[i]->bcDirection(); - secondBoundaryCache_[i] = points_[i]->secondBoundaryConds(); - corrCache_[i] = points_[i]->corrPoint(); + bcTypeCache_[i] = points_[i].bcType(); + directionCache_[i] = points_[i].bcDirection(); + secondBoundaryCache_[i] = points_[i].secondBoundaryConds(); + corrCache_[i] = points_[i].corrPoint(); } totalCoefMatrix_.resize(nTotalPoints_, nTotalPoints_); @@ -175,13 +175,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:{ @@ -232,16 +232,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); @@ -289,11 +289,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:{ @@ -309,31 +309,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; @@ -366,8 +366,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: { @@ -379,8 +379,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: { @@ -394,8 +394,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]; @@ -407,8 +407,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: { @@ -425,8 +425,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: { @@ -437,8 +437,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: { @@ -462,54 +462,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. @@ -551,8 +551,8 @@ namespace FVM_ANDS{ isSouthBoundary = direction == FaceDirection::SOUTH; //Corner cases... - bool secondaryWestBound = (secondBoundaryCache_[i]->direction == FaceDirection::WEST); - bool secondaryEastBound = (secondBoundaryCache_[i]->direction == FaceDirection::EAST); + bool secondaryWestBound = (secondBoundaryCache_[i].direction == FaceDirection::WEST); + bool secondaryEastBound = (secondBoundaryCache_[i].direction == FaceDirection::EAST); isWestBoundary = (direction == FaceDirection::WEST || secondaryWestBound); isEastBoundary = (direction == FaceDirection::EAST || secondaryEastBound); @@ -560,8 +560,8 @@ namespace FVM_ANDS{ //only call this lookup function on boundary nodes which are inconsequential in number idx_N = isNorthBoundary? corrCache_[i] : idx_N; idx_S = isSouthBoundary? corrCache_[i] : idx_S; - idx_E = isEastBoundary? (secondaryEastBound ? secondBoundaryCache_[i]->corrPoint : corrCache_[i]) : idx_E; - idx_W = isWestBoundary? (secondaryEastBound ? secondBoundaryCache_[i]->corrPoint : corrCache_[i]) : idx_W; + idx_E = isEastBoundary? (secondaryEastBound ? secondBoundaryCache_[i].corrPoint : corrCache_[i]) : idx_E; + idx_W = isWestBoundary? (secondaryEastBound ? secondBoundaryCache_[i].corrPoint : corrCache_[i]) : idx_W; } //When you declare these vars (inside or outside loop) has 0 impact) //takes ~ 6 out of 18 ns on background var calcs @@ -592,7 +592,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]); @@ -601,7 +601,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]); @@ -611,7 +611,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]); @@ -621,7 +621,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]); From a7278c29a91067315ec81006d874507abcd50293 Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Tue, 30 Dec 2025 16:51:03 +0800 Subject: [PATCH 5/8] fix build errors --- Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp | 6 +++--- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 11 +++-------- 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp index 2b3d4cc0..9ea55e95 100644 --- a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp +++ b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp @@ -6,9 +6,9 @@ #include "FVM_ANDS_HelperFunctions.hpp" #include #include -using PointVariant = std::variant; 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); @@ -150,7 +150,7 @@ namespace FVM_ANDS{ std::vector points_; std::vector bcCache_; // Cache bcType() std::vector directionCache_; // Cache bcDirection() - std::vector> secondBounaryCache_; // Cache secondBoundaryConds() + std::vector> secondBoundaryCache_; // Cache secondBoundaryConds() std::vector corrCache_; // Cache corrPoint() Eigen::SparseMatrix totalCoefMatrix_; Eigen::VectorXd rhs_; @@ -323,7 +323,7 @@ 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(); + Point* point = points_[pointID]; if(point->bcType() == BoundaryConditionFlag::INTERIOR) return neighbor_point_interior(direction, pointID); if(point->bcDirection() == direction){ diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index 058012a6..5182ff1d 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) : @@ -49,7 +50,7 @@ namespace FVM_ANDS{ return Point(); // Construct directly in vector }); for (size_t i = 0; i < points_.size(); ++i) { - bcTypeCache_[i] = points_[i].bcType(); + bcCache_[i] = points_[i].bcType(); directionCache_[i] = points_[i].bcDirection(); secondBoundaryCache_[i] = points_[i].secondBoundaryConds(); corrCache_[i] = points_[i].corrPoint(); @@ -517,12 +518,6 @@ namespace FVM_ANDS{ Eigen::VectorXd AdvDiffSystem::forwardEulerAdvection(bool operatorSplit, bool parallelAdvection) const noexcept{ Eigen::VectorXd soln(nTotalPoints_); - // cache results outside of hotloop - for (const auto& point : points_) { - bcTypes_.push_back(std::visit([](auto&& p) { return p.bcType(); }, point)); - directions_.push_back(std::visit([](auto&& p) { return p.bcDirection(); }, point)); - } - // double avgBackgroundCalcTime = 0; //Explicit Time-Stepping #pragma omp parallel for \ @@ -546,7 +541,7 @@ namespace FVM_ANDS{ // isNorthBoundary = direction == FaceDirection::NORTH; // isSouthBoundary = direction == FaceDirection::SOUTH; - FaceDirection direction = directions_[i]; + FaceDirection direction = directionCache_[i]; isNorthBoundary = direction == FaceDirection::NORTH; isSouthBoundary = direction == FaceDirection::SOUTH; From 2e4214d02545a3147c3dffe6a502bdd41c867df3 Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Tue, 30 Dec 2025 16:57:08 +0800 Subject: [PATCH 6/8] fixing build errors --- Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp | 10 +++++----- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 10 ++++++---- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp index 9ea55e95..d2ea5b8b 100644 --- a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp +++ b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp @@ -323,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]; - 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){ + 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/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index 5182ff1d..6e878e86 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -546,8 +546,10 @@ namespace FVM_ANDS{ isSouthBoundary = direction == FaceDirection::SOUTH; //Corner cases... - bool secondaryWestBound = (secondBoundaryCache_[i].direction == FaceDirection::WEST); - bool secondaryEastBound = (secondBoundaryCache_[i].direction == FaceDirection::EAST); + if (secondBoundaryCache_[i]) { + bool secondaryWestBound = (secondBoundaryCache_[i]->direction == FaceDirection::WEST); + bool secondaryEastBound = (secondBoundaryCache_[i]->direction == FaceDirection::EAST); + } isWestBoundary = (direction == FaceDirection::WEST || secondaryWestBound); isEastBoundary = (direction == FaceDirection::EAST || secondaryEastBound); @@ -555,8 +557,8 @@ namespace FVM_ANDS{ //only call this lookup function on boundary nodes which are inconsequential in number idx_N = isNorthBoundary? corrCache_[i] : idx_N; idx_S = isSouthBoundary? corrCache_[i] : idx_S; - idx_E = isEastBoundary? (secondaryEastBound ? secondBoundaryCache_[i].corrPoint : corrCache_[i]) : idx_E; - idx_W = isWestBoundary? (secondaryEastBound ? secondBoundaryCache_[i].corrPoint : corrCache_[i]) : idx_W; + idx_E = isEastBoundary? (secondaryEastBound ? secondBoundaryCache_[i]->corrPoint : corrCache_[i]) : idx_E; + idx_W = isWestBoundary? (secondaryEastBound ? secondBoundaryCache_[i]->corrPoint : corrCache_[i]) : idx_W; } //When you declare these vars (inside or outside loop) has 0 impact) //takes ~ 6 out of 18 ns on background var calcs From d92d964f327e56179ec9c6e6449397ab5bd23e5a Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Tue, 30 Dec 2025 16:58:16 +0800 Subject: [PATCH 7/8] fixing build errors --- Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp index d2ea5b8b..050928cb 100644 --- a/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp +++ b/Code.v05-00/include/FVM_ANDS/AdvDiffSystem.hpp @@ -327,7 +327,7 @@ namespace FVM_ANDS{ if(point.bcType() == BoundaryConditionFlag::INTERIOR) return neighbor_point_interior(direction, pointID); if(point.bcDirection() == direction){ - return point->corrPoint(); + return point.corrPoint(); } else if (point.secondBoundaryConds() && point.secondBoundaryConds().value().direction == direction){ return point.secondBoundaryConds().value().corrPoint; From 3d7274888aade4e7d5ca2f83211a29b532167d90 Mon Sep 17 00:00:00 2001 From: Coco Yeung Date: Tue, 30 Dec 2025 10:04:57 +0000 Subject: [PATCH 8/8] Fix all build errors --- Code.v05-00/include/Core/LAGRIDPlumeModel.hpp | 4 +-- .../include/FVM_ANDS/BoundaryCondition.hpp | 4 +-- Code.v05-00/include/FVM_ANDS/FVM_Solver.hpp | 2 +- Code.v05-00/src/Core/LAGRIDPlumeModel.cpp | 2 +- Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp | 26 ++++++++++++++----- Code.v05-00/tests/test_adv_diff_solver.cpp | 6 ++--- 6 files changed, 29 insertions(+), 15 deletions(-) diff --git a/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp b/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp index b4968c42..4060cd08 100644 --- a/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp +++ b/Code.v05-00/include/Core/LAGRIDPlumeModel.hpp @@ -114,8 +114,8 @@ class LAGRIDPlumeModel { const std::vector>& mask, Vector_1D& xEdgesNew, Vector_1D& yEdgesNew, - Vector_1D& xCoordsNew, V - ector_1D& yCoordsNew + Vector_1D& xCoordsNew, + Vector_1D& yCoordsNew ); Vector_2D applyWeights( const Eigen::SparseMatrix& weights, diff --git a/Code.v05-00/include/FVM_ANDS/BoundaryCondition.hpp b/Code.v05-00/include/FVM_ANDS/BoundaryCondition.hpp index fc6d0b85..8027c055 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 27c78ca6..9da52927 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/Core/LAGRIDPlumeModel.cpp b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp index eb36eca1..bcb0db4e 100644 --- a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp +++ b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp @@ -367,7 +367,7 @@ void LAGRIDPlumeModel::updateDiffVecs() { } } } -void LAGRIDPlumeModel::runTrans,port(double timestep) { +void LAGRIDPlumeModel::runTransport(double timestep) { //Update the zero bc to reflect grid size changes auto ZERO_BC = FVM_ANDS::bcFrom2DVector(iceAerosol_.getPDF()[0], true); diff --git a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp index 6e878e86..861fede9 100644 --- a/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp +++ b/Code.v05-00/src/FVM_ANDS/AdvDiffSystem.cpp @@ -544,11 +544,16 @@ namespace FVM_ANDS{ 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... - if (secondBoundaryCache_[i]) { - bool secondaryWestBound = (secondBoundaryCache_[i]->direction == FaceDirection::WEST); - bool secondaryEastBound = (secondBoundaryCache_[i]->direction == FaceDirection::EAST); + if (secondBC) { + secondaryWestBound = (secondBC->direction == FaceDirection::WEST); + secondaryEastBound = (secondBC->direction == FaceDirection::EAST); } isWestBoundary = (direction == FaceDirection::WEST || secondaryWestBound); @@ -557,8 +562,17 @@ namespace FVM_ANDS{ //only call this lookup function on boundary nodes which are inconsequential in number idx_N = isNorthBoundary? corrCache_[i] : idx_N; idx_S = isSouthBoundary? corrCache_[i] : idx_S; - idx_E = isEastBoundary? (secondaryEastBound ? secondBoundaryCache_[i]->corrPoint : corrCache_[i]) : idx_E; - idx_W = isWestBoundary? (secondaryEastBound ? secondBoundaryCache_[i]->corrPoint : corrCache_[i]) : idx_W; + 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 diff --git a/Code.v05-00/tests/test_adv_diff_solver.cpp b/Code.v05-00/tests/test_adv_diff_solver.cpp index 286fdee5..219b6cbc 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);