From 7a1d9bb5b7ac4e01b5bfddcf748d16232dcf6be6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 5 Feb 2025 15:42:36 +0100 Subject: [PATCH] feat: Implement steady-state hydraulic charge calculation --- include/tug/Advection/Velocities.hpp | 160 +++++++++++++------ include/tug/Core/Numeric/FTCS.hpp | 20 ++- include/tug/Core/Numeric/SimulationInput.hpp | 3 +- include/tug/tug.hpp | 1 - test/testDiffusion.cpp | 33 +--- test/testVelocities.cpp | 24 +-- 6 files changed, 133 insertions(+), 108 deletions(-) diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp index dd495fa..2fbd774 100644 --- a/include/tug/Advection/Velocities.hpp +++ b/include/tug/Advection/Velocities.hpp @@ -8,6 +8,7 @@ #pragma once +#include "tug/Core/Numeric/SimulationInput.hpp" #include #include #include @@ -16,12 +17,12 @@ #include #include +#include #include #include #include #include #include -#include #ifdef _OPENMP #include @@ -32,38 +33,38 @@ using namespace Eigen; namespace tug { -enum HYDRAULIC_MODE { TRANSIENT, STEADY_STATE }; -enum HYDRAULIC_RESOLVE { EXPLICIT, IMPLICIT }; +enum class HYDRAULIC_MODE { TRANSIENT, STEADY_STATE }; +enum class HYDRAULIC_RESOLVE { EXPLICIT, IMPLICIT }; template -class Velocities : public BaseSimulation { +class Velocities : public BaseSimulationGrid { private: int innerIterations{1}; T timestep{-1}; T epsilon{1E-5}; int numThreads{omp_get_num_procs()}; - Grid &grid; - Boundary &bc; - RowMajMat velocitiesX; RowMajMat velocitiesY; + RowMajMat alphaX; + RowMajMat alphaY; + public: - /** - * @brief Construct a new Velocities object, used to calculate Hydraulic - * Charge and Darcy-Velocities. A timestep and a number of iterations can be - * set. By default iterations is set to -1. If the number of iterations is set - * to a value below 1 the simulation will run until the Hydraulic Charge - * converges. The Epsilon value to check convergence can be set, the default - * is 1E-5. CSV Output is off by default. - * - * @param grid Valid grid object - * @param bc Valid boundary condition object - */ - Velocities(Grid &_grid, Boundary &_bc) - : grid(_grid), bc(_bc), velocitiesX(grid.getRow(), grid.getCol() + 1), - velocitiesY(grid.getRow() + 1, grid.getCol()) {}; + Velocities(RowMajMat &origin) + : BaseSimulationGrid(origin), + velocitiesX(origin.rows(), origin.cols() + 1), + velocitiesY(origin.rows() + 1, origin.cols()), + alphaX(origin.rows(), origin.cols()), + alphaY(origin.rows(), origin.cols()) {}; + + Velocities(T *data, std::size_t rows, std::size_t cols) + : BaseSimulationGrid(data, rows, cols), velocitiesX(rows, cols + 1), + velocitiesY(rows + 1, cols), alphaX(rows, cols), alphaY(rows, cols) {}; + + Velocities(T *data, std::size_t length) + : BaseSimulationGrid(data, 1, length), velocitiesX(1, length + 1), + alphaX(1, length) {}; /** * @brief Set the epsilon value, the relativ error allowed for convergence @@ -77,22 +78,61 @@ public: this->epsilon = epsilon; } + /** + * @brief Get the alphaX matrix. + * + * @return RowMajMat& Reference to the alphaX matrix. + */ + RowMajMat &getAlphaX() { return alphaX; } + + /** + * @brief Get the alphaY matrix. + * + * @return RowMajMat& Reference to the alphaY matrix. + */ + RowMajMat &getAlphaY() { + tug_assert( + this->getDim(), + "Grid is not two dimensional, there is no domain in y-direction!"); + + return alphaY; + } + + /** + * @brief Set the alphaX matrix. + * + * @param alphaX The new alphaX matrix. + */ + void setAlphaX(const RowMajMat &alphaX) { this->alphaX = alphaX; } + + /** + * @brief Set the alphaY matrix. + * + * @param alphaY The new alphaY matrix. + */ + void setAlphaY(const RowMajMat &alphaY) { + tug_assert( + this->getDim(), + "Grid is not two dimensional, there is no domain in y-direction!"); + + this->alphaY = alphaY; + } + /** * @brief Set the timestep per iteration * * @param timestep timestep per iteration */ - void setTimestep(T timestep) { + void setTimestep(T timestep) override { if (timestep <= 0) { throw std::invalid_argument("Timestep must be greater than zero"); } this->timestep = timestep; - const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + const T deltaColSquare = this->deltaCol() * this->deltaCol(); + const T deltaRowSquare = this->deltaRow() * this->deltaRow(); const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); - const T maxK = - std::max(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff()); + const T maxK = std::max(this->alphaX.maxCoeff(), this->alphaY.maxCoeff()); T cfl = minDeltaSquare / (4 * maxK); if (timestep > cfl) { @@ -146,32 +186,45 @@ public: * @brief Simulation of hydraulic charge either until convergence, * or for a number of selected timesteps. Calculation of Darcy-velocities. */ - void run() { + void run() override { // if iterations < 1 calculate hydraulic charge until steady state is // reached - if constexpr (hyd_mode == STEADY_STATE) { - // Calculate largest possible timestep, depending on K and gridsize - const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + SimulationInput input = {.concentrations = + this->getConcentrationMatrix(), + .alphaX = this->getAlphaX(), + .alphaY = this->getAlphaY(), + .boundaries = this->getBoundaryConditions(), + .dim = this->getDim(), + .timestep = this->timestep, + .rowMax = this->rows(), + .colMax = this->cols(), + .deltaRow = this->deltaRow(), + .deltaCol = this->deltaCol()}; + + if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) { + const T deltaColSquare = this->deltaCol() * this->deltaCol(); + const T deltaRowSquare = this->deltaRow() * this->deltaRow(); const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); - const T maxK = - std::max(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff()); + const T maxK = std::max(this->alphaX.maxCoeff(), this->alphaY.maxCoeff()); + // Calculate largest possible timestep, depending on K and gridsize setTimestep(minDeltaSquare / (4 * maxK)); + input.timestep = this->timestep; + RowMajMat oldConcentrations; do { - oldConcentrations = grid.getConcentrations(); - (void)calculate_hydraulic_flow(); - } while (!checkConvergance(oldConcentrations, grid.getConcentrations())); + oldConcentrations = this->getConcentrationMatrix(); + (void)calculate_hydraulic_flow(input); + } while (!checkConvergance(oldConcentrations)); } else { if (timestep == -1) { throw_invalid_argument("Timestep is not set"); } for (int i = 0; i < innerIterations; i++) { - (void)calculate_hydraulic_flow(); + (void)calculate_hydraulic_flow(input); } } @@ -182,11 +235,11 @@ private: /** * @brief Calculate the new hydraulic charge using FTCS */ - void calculate_hydraulic_flow() { - if constexpr (hyd_resolve == EXPLICIT) { - FTCS_2D(this->grid, this->bc, this->timestep, this->numThreads); + void calculate_hydraulic_flow(SimulationInput &sim_in) { + if constexpr (hyd_resolve == HYDRAULIC_RESOLVE::EXPLICIT) { + FTCS_2D(sim_in, numThreads); } else { - BTCS_2D(this->grid, this->bc, this->timestep, ThomasAlgorithm); + BTCS_2D(sim_in, ThomasAlgorithm, numThreads); } }; @@ -197,10 +250,9 @@ private: * containing old and new Charge values, the relative error is below the * selected Epsilon */ - bool checkConvergance(const RowMajMat &oldHeads, - const RowMajMat &newHeads) { - const auto abs_err = (oldHeads - newHeads).cwiseAbs(); - const auto rel_err = abs_err.cwiseQuotient(newHeads); + bool checkConvergance(const RowMajMat &oldHeads) { + const auto abs_err = (oldHeads - this->getConcentrationMatrix()).cwiseAbs(); + const auto rel_err = abs_err.cwiseQuotient(this->getConcentrationMatrix()); return rel_err.maxCoeff() < epsilon; } @@ -210,14 +262,16 @@ private: * directions */ void computeFluidVelocities() { - const std::size_t rows = grid.getRow(); - const std::size_t cols = grid.getCol(); - const T dx = grid.getDeltaRow(); - const T dy = grid.getDeltaCol(); - const RowMajMat &hydraulicCharges = grid.getConcentrations(); + const std::size_t rows = this->rows(); + const std::size_t cols = this->cols(); + const T dx = this->deltaCol(); + const T dy = this->deltaRow(); + const RowMajMat &hydraulicCharges = this->getConcentrationMatrix(); - const RowMajMat &permKX = grid.getAlphaX(); - const RowMajMat &permKY = grid.getAlphaY(); + const RowMajMat &permKX = this->alphaX; + const RowMajMat &permKY = this->alphaY; + + const Boundary &bc = this->getBoundaryConditions(); // calculate velocities in x-direction for (std::size_t i_rows = 0; i_rows < rows; i_rows++) { @@ -229,7 +283,7 @@ private: } case BC_TYPE_CONSTANT: { velocitiesX(i_rows, 0) = - -permKX(i_rows, 0) * + -this->alphaX(i_rows, 0) * (hydraulicCharges(i_rows, 0) - bc_left.getValue()) / (dx / 2); break; } diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp index 3ab5cc2..bcb78ba 100644 --- a/include/tug/Core/Numeric/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -8,13 +8,13 @@ #ifndef FTCS_H_ #define FTCS_H_ -#include "tug/Core/TugUtils.hpp" +#include "tug/Core/Matrix.hpp" #include #include -#include #include #include #include +#include #ifdef _OPENMP #include @@ -57,17 +57,15 @@ constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center, } template -static inline void checkAndSetConstantInnerCells(const Boundary &bc, - Grid &grid) { +static inline void +checkAndSetConstantInnerCells(const Boundary &bc, + RowMajMatMap &concentrations, std::size_t rows, + std::size_t cols) { const auto &inner_boundaries = bc.getInnerBoundaries(); if (inner_boundaries.empty()) { return; } - auto &concentrations = grid.getConcentrations(); - const auto rows = grid.getRow(); - const auto cols = grid.getCol(); - for (const auto &[rowcol, value] : inner_boundaries) { const auto &row = rowcol.first; const auto &col = rowcol.second; @@ -90,6 +88,9 @@ template static void FTCS_1D(SimulationInput &input) { const auto &alphaX = input.alphaX; const auto &bc = input.boundaries; + checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax, + input.colMax); + // only one row in 1D case -> row constant at index 0 int row = 0; @@ -168,6 +169,9 @@ static void FTCS_2D(SimulationInput &input, int numThreads) { const auto &alphaY = input.alphaY; const auto &bc = input.boundaries; + checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax, + input.colMax); + const T sx = timestep / (deltaCol * deltaCol); const T sy = timestep / (deltaRow * deltaRow); diff --git a/include/tug/Core/Numeric/SimulationInput.hpp b/include/tug/Core/Numeric/SimulationInput.hpp index d68effc..5068537 100644 --- a/include/tug/Core/Numeric/SimulationInput.hpp +++ b/include/tug/Core/Numeric/SimulationInput.hpp @@ -11,9 +11,8 @@ template struct SimulationInput { const RowMajMat &alphaX; const RowMajMat &alphaY; const Boundary boundaries; - const std::uint8_t dim; - const T timestep; + T timestep; const std::size_t rowMax; const std::size_t colMax; const T deltaRow; diff --git a/include/tug/tug.hpp b/include/tug/tug.hpp index 89440b9..bc797ae 100644 --- a/include/tug/tug.hpp +++ b/include/tug/tug.hpp @@ -4,4 +4,3 @@ #include #include #include -#include \ No newline at end of file diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 9a74670..75a51b7 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -3,8 +3,7 @@ #include "gtest/gtest.h" #include #include - -#include +#include #include #include @@ -215,33 +214,3 @@ DIFFUSION_TEST(ConstantInnerCell) { EXPECT_FALSE((concentrations_result.array() < 0.0).any()); } - -DIFFUSION_TEST(ConstantInnerCellFTCS) { - constexpr std::uint32_t nrows = 5; - constexpr std::uint32_t ncols = 5; - - auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0); - auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - - tug::Grid64 grid(concentrations); - grid.setAlpha(alphax, alphay); - - tug::Boundary bc(grid); - // inner - bc.setInnerBoundary(2, 2, 0); - - tug::Diffusion sim(grid, bc); - sim.setTimestep(1); - sim.setIterations(1); - - MatrixXd input_values(concentrations); - sim.run(); - - EXPECT_DOUBLE_EQ(grid.getConcentrations()(2, 2), 0); - EXPECT_LT(grid.getConcentrations().sum(), input_values.sum()); - - EXPECT_FALSE((grid.getConcentrations().array() > 1.0).any()); - - EXPECT_FALSE((grid.getConcentrations().array() < 0.0).any()); -} \ No newline at end of file diff --git a/test/testVelocities.cpp b/test/testVelocities.cpp index 03b1450..9e4ecac 100644 --- a/test/testVelocities.cpp +++ b/test/testVelocities.cpp @@ -25,11 +25,15 @@ VELOCITIES_TEST(SteadyStateCenter) { tug::RowMajMat permKY = tug::RowMajMat::Constant(rows, cols, K); - tug::Grid64 gridHeads(hydHeads); - gridHeads.setDomain(100, 100); - gridHeads.setAlpha(permKX, permKY); + tug::Velocities + velo(hydHeads); - tug::Boundary bcH = tug::Boundary(gridHeads); + velo.setDomain(100, 100); + velo.setAlphaX(permKX); + velo.setAlphaY(permKY); + + tug::Boundary &bcH = velo.getBoundaryConditions(); bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1); bcH.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1); bcH.setBoundarySideConstant(tug::BC_SIDE_TOP, 1); @@ -37,14 +41,10 @@ VELOCITIES_TEST(SteadyStateCenter) { bcH.setInnerBoundary(center_row, center_col, 10); - tug::Velocities - velocities(gridHeads, bcH); + velo.run(); - velocities.run(); - - const auto &velocitiesX = velocities.getVelocitiesX(); - const auto &velocitiesY = velocities.getVelocitiesY(); + const auto &velocitiesX = velo.getVelocitiesX(); + const auto &velocitiesY = velo.getVelocitiesY(); // Assert @@ -81,4 +81,4 @@ VELOCITIES_TEST(SteadyStateCenter) { } } } -} \ No newline at end of file +}