From 763a17b80f1bb464b489baa86f087353e36daa0d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 10 Dec 2024 11:52:03 +0100 Subject: [PATCH 01/20] feat: Add implementation of Advection from Christopher Eschenbach (does not work!) --- examples/Advection.cpp | 59 ++++ examples/CMakeLists.txt | 2 + include/tug/Advection.hpp | 368 ++++++++++++++++++++++++ include/tug/Core/Velocities.hpp | 488 ++++++++++++++++++++++++++++++++ include/tug/Matrix.hpp | 0 5 files changed, 917 insertions(+) create mode 100644 examples/Advection.cpp create mode 100644 include/tug/Advection.hpp create mode 100644 include/tug/Core/Velocities.hpp create mode 100644 include/tug/Matrix.hpp diff --git a/examples/Advection.cpp b/examples/Advection.cpp new file mode 100644 index 0000000..414a8af --- /dev/null +++ b/examples/Advection.cpp @@ -0,0 +1,59 @@ +#include +#include +#include +#include + +using namespace Eigen; +using namespace tug; + +int main(int argc, char *argv[]) { + int row = 21; + int col = 21; + // create two grids of equal size, grid1 for hydraulics heads, grid2 for + // Concentrations + + RowMajMat hydHeads = RowMajMat::Constant(row, col, 1); + RowMajMat concentrations = RowMajMat::Constant(row, col, 0); + + Grid64 gridHeads(hydHeads); + Grid64 gridConc(concentrations); + gridHeads.setDomain(100, 100); + gridConc.setDomain(100, 100); + + // create boundaries + Boundary bcH = Boundary(gridHeads); + bcH.setBoundarySideConstant(BC_SIDE_LEFT, 10); + bcH.setBoundarySideConstant(BC_SIDE_RIGHT, 1); + bcH.setBoundarySideClosed(BC_SIDE_TOP); + bcH.setBoundarySideClosed(BC_SIDE_BOTTOM); + Boundary bcC = Boundary(gridConc); + bcC.setBoundarySideConstant(BC_SIDE_LEFT, 0.1); + bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 1); + bcC.setBoundarySideClosed(BC_SIDE_TOP); + bcC.setBoundarySideClosed(BC_SIDE_BOTTOM); + + Velocities velocities = Velocities(gridHeads, bcH); + velocities.setOutputCSV(CSV_OUTPUT_ON); + velocities.setK(1E-2); + velocities.setEpsilon(1E-8); + velocities.setInjh(10); + velocities.setIterations(0); + // calculate steady hydraulic heads + velocities.run(); + + std::cout << "Velocities simulation finished." << std::endl; + std::cout << hydHeads << std::endl; + + // set true for steady case + Advection advection = Advection(velocities, gridConc, bcC, true); + advection.setPorosity(0.2); + advection.setIterations(21); + // set timestep close almost exactly to clf to test advection + advection.setTimestep(5039.05); + // velocities.setOutputCSV(CSV_OUTPUT_VERBOSE); + advection.run(); + + std::cout << "Advection simulation finished." << std::endl; + + std::cout << concentrations << std::endl; +} \ No newline at end of file diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 3c8bf28..8924b9c 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,7 +1,9 @@ add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp) add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp) add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp) +add_executable(Advection Advection.cpp) target_link_libraries(BTCS_2D_proto_example tug) target_link_libraries(FTCS_2D_proto_closed_mdl tug) target_link_libraries(FTCS_2D_proto_example_mdl tug) +target_link_libraries(Advection tug) diff --git a/include/tug/Advection.hpp b/include/tug/Advection.hpp new file mode 100644 index 0000000..4d530bf --- /dev/null +++ b/include/tug/Advection.hpp @@ -0,0 +1,368 @@ +/** + * @file Advection.hpp + * @brief API of Advection class, holding information for a simulation of + * advection. Holds a predifined Grid object, Boundary object and Velocities + * object + */ + +#pragma once + +#include "tug/Core/Matrix.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +using namespace Eigen; +namespace tug { +template class Advection { +public: + /** + * @brief Construct a new Advection object, used to calculate material + * transport. A timestep and number of iterations must be set. A transient + * case can be selected by initializing Steady=false. With each timestep the + * Velocities object will also be updated. + * A steady case can be selected by initializing Steady=true. The velocities + * object will not be updated. Velocities can be simulated to convergence + * beforehand. Porosity can be set, the default is 1. CSV Output is off by + * default. + * + * @param grid Valid grid object + * @param bc Valid Boundary object + * @param Steady Used to choose between Steady and Transient case. Either true + * or false + */ + Advection(Velocities &_velocities, Grid &_grid, Boundary &_bc, + bool Steady) + : velocities(_velocities), grid(_grid), bc(_bc), + outx(_velocities.getOutx()), outy(_velocities.getOuty()), + Steady(Steady) {}; + + /** + * @brief Sets the porosity of the medium + * + * @param porosity new porosity value + */ + void setPorosity(T porosity) { + if (porosity < 0 || porosity > 1) { + throw std::invalid_argument( + "Porosity must be a value between 0 and 1 (inclusive)"); + } + this->porosity = porosity; + } + + /** + * @brief Set the desired iterations to be calculated. A value greater + * than zero must be specified here. Setting iterations is required. + * + * @param iterations Number of iterations to be simulated. + */ + void setIterations(int iterations) { + if (iterations <= 0) { + throw std::invalid_argument( + "Number of iterations must be greater than zero. Provided value: " + + std::to_string(iterations)); + } + this->iterations = iterations; + }; + + /** + * @brief Set the size of the timestep. Must be greater than zero + * + * @param timestep Size of the timestep + */ + void setTimestep(T timestep) { + if (timestep <= 0) { + throw std::invalid_argument( + "Timestep must be greater than zero. Provided value: " + + std::to_string(timestep)); + } else { + this->timestep = timestep; + } + } + + /** + * @brief Set the number of desired openMP Threads. + * + * @param num_threads Number of desired threads. Must have a value between + * 1 and the maximum available number of processors. The + * maximum number of processors is set as the default case during Advection + * construction. + */ + void setNumberThreads(int num_threads) { + if (num_threads > 0 && num_threads <= omp_get_num_procs()) { + this->numThreads = num_threads; + } else { + int maxThreadNumber = omp_get_num_procs(); + if (num_threads > maxThreadNumber) { + throw std::invalid_argument( + "Number of threads exceeds the number of processor cores (" + + std::to_string(maxThreadNumber) + ")."); + } else { + throw std::invalid_argument("Number of threads is less than 1."); + } + } + } + + /** + * @brief Set the option to output the results to a CSV file. Off by default. + * + * + * @param csv_output Valid output option. The following options can be set + * here: + * - CSV_OUTPUT_OFF: do not produce csv output + * - CSV_OUTPUT_ON: produce csv output with last + * concentration matrix + if (csv_output < CSV_OUTPUT_OFF || csv_output > CSV_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid CSV output option given: " + + std::to_string(csv_output)); + } + void setOutputCSV(CSV_OUTPUT csv_output) { + if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid CSV output option given!"); + } + + this->csv_output = csv_output; + } + + /** + * @brief Creates a CSV file with a name containing the current simulation + * parameters. If the data name already exists, an additional counter + * is appended to the name. The name of the file is built up as follows: + * + + + + * +.csv + * + * @return string Filename with configured simulation parameters. + */ + std::string createCSVfile(std::string Type) const { + std::ofstream file; + int appendIdent = 0; + std::string appendIdentString; + + std::string row = std::to_string(grid.getRow()); + std::string col = std::to_string(grid.getCol()); + std::string numIterations = std::to_string(iterations); + + std::string filename = + Type + "_" + row + "_" + col + "_" + numIterations + ".csv"; + + while (std::filesystem::exists(filename)) { + appendIdent += 1; + appendIdentString = std::to_string(appendIdent); + filename = Type + "_" + row + "_" + col + "_" + numIterations + "-" + + appendIdentString + ".csv"; + } + + file.open(filename); + if (!file) { + throw std::runtime_error("Failed to open file: " + filename); + } + + file.close(); + + return filename; + } + /** + * @brief Writes the currently calculated Concentration values of the grid + * into the CSV file with the passed filename. + * + * @param filename Name of the file to which the Concentration values are + * to be written. + */ + void printConcentrationCSV(const std::string &filename) const { + std::ofstream file; + + file.open(filename, std::ios_base::app); + if (!file) { + throw std::runtime_error("Failed to open file: " + filename); + } + + Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); + file << grid.getConcentrations().format(do_not_align) << std::endl; + file << std::endl << std::endl; + file.close(); + } + + /** + * @brief function calculating material transport for one timestep + */ + void adv() { + int rows = grid.getRow(); + int cols = grid.getCol(); + T volume = grid.getDeltaRow() * grid.getDeltaCol(); + + RowMajMat &newConcentrations = grid.getConcentrations(); + + // Calculate Courant-Levy-Frederich condition + T maxFx = std::max(abs(outx.maxCoeff()), abs(outx.minCoeff())); + T maxFy = std::max(abs(outy.maxCoeff()), abs(outy.minCoeff())); + T maxF = std::max(maxFx, maxFy); + + if (maxF == 0) { + throw std::runtime_error("Division by zero: maxF is zero."); + } + + T cvf = abs((volume * porosity) / maxF); + int innerSteps = (int)ceil(timestep / cvf); + T innerTimestep = timestep / innerSteps; + + for (int k = 0; k < innerSteps; k++) { + const Eigen::MatrixX oldConcentrations = newConcentrations; +// Calculate sum of incoming/outgoing Flow*Concentration in x-direction in each +// cell +#pragma omp parallel for num_threads(numThreads) schedule(static) + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols + 1; j++) { + if (j == 0) { + if (bc.getBoundaryElementType(BC_SIDE_LEFT, i) != BC_TYPE_CLOSED) { + if (outx(i, j) > 0) { + // outx positive -> flow from border to cell i,j + newConcentrations(i, j) += + outx(i, j) * bc.getBoundaryElementValue(BC_SIDE_LEFT, i); + } else if (outx(i, j) < 0) { + // outx negative -> flow from i,j towards border + newConcentrations(i, j) += outx(i, j) * oldConcentrations(i, j); + } + } + } else if (j == cols) { + if (bc.getBoundaryElementType(BC_SIDE_RIGHT, i) != BC_TYPE_CLOSED) { + if (outx(i, j) > 0) { + // outx positive-> flow from i,j-1 towards border + newConcentrations(i, j - 1) -= + outx(i, j) * oldConcentrations(i, j - 1); + } else if (outx(i, j) < 0) { + // outx negative -> flow from border to cell i,j-1 + newConcentrations(i, j - 1) -= + outx(i, j) * bc.getBoundaryElementValue(BC_SIDE_LEFT, i); + } + } + } + // flow between inner cells + else { + // outx positive -> flow from cell i,j-1 towards cell i,j + if (outx(i, j) > 0) { + newConcentrations(i, j - 1) -= + outx(i, j) * oldConcentrations(i, j - 1); + newConcentrations(i, j) += + outx(i, j) * oldConcentrations(i, j - 1); + } + // outx negative -> flow from cell i,j toward cell i,j-1 + else if (outx(i, j) < 0) { + newConcentrations(i, j - 1) -= + outx(i, j) * oldConcentrations(i, j); + newConcentrations(i, j) += outx(i, j) * oldConcentrations(i, j); + } + } + } + } +// calculate sum in y-direction +// parallelize outer loop over columns to ensure thread-safety, each thread only +// modifies cells within its column +#pragma omp parallel for num_threads(numThreads) + for (int j = 0; j < cols; j++) { + for (int i = 0; i < rows + 1; i++) { + if (i == 0) { + if (bc.getBoundaryElementType(BC_SIDE_TOP, j) != BC_TYPE_CLOSED) { + if (outy(i, j) > 0) { + // outy positive -> flow from border to cell i,j + newConcentrations(i, j) += + outy(i, j) * bc.getBoundaryElementValue(BC_SIDE_TOP, j); + } else if (outy(i, j) < 0) { + // outy negative -> flow from i,j towards border + newConcentrations(i, j) += outy(i, j) * oldConcentrations(i, j); + } + } + } else if (i == rows) { + if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, j) != + BC_TYPE_CLOSED) { + if (outy(i, j) > 0) { + // outy positive-> flow from i-1,j towards border + newConcentrations(i - 1, j) -= + outy(i, j) * oldConcentrations(i - 1, j); + } else if (outy(i, j) < 0) { + // outy negative -> flow from border to cell i,j-1 + newConcentrations(i - 1, j) -= + outy(i, j) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, j); + } + } + } + // flow between inner cells + else { + // outy positive -> flow from cell i-1,j towards cell i,j + if (outy(i, j) > 0) { + newConcentrations(i - 1, j) -= + outy(i, j) * oldConcentrations(i - 1, j); + newConcentrations(i, j) += + outy(i, j) * oldConcentrations(i - 1, j); + } + // outy negative -> flow from cell i,j toward cell i-1,j + else if (outy(i, j) < 0) { + newConcentrations(i - 1, j) -= + outy(i, j) * oldConcentrations(i, j); + newConcentrations(i, j) += outy(i, j) * oldConcentrations(i, j); + } + } + } + } + newConcentrations = + oldConcentrations + + newConcentrations * (innerTimestep / (porosity * volume)); + } + } + + void run() { + std::string filename; + if (csv_output >= CSV_OUTPUT_ON) { + filename = createCSVfile("Concentrations"); + } + + if (Steady == false) { + velocities.setTimestep(timestep); + velocities.setIterations(1); + } + + for (int i = 0; i < iterations; i++) { + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printConcentrationCSV(filename); + } + // if steady==false update charge and velocities with equal timestep + if (Steady == false) { + velocities.run(); + } + adv(); + } + + if (csv_output >= CSV_OUTPUT_ON) { + printConcentrationCSV(filename); + } + } + +private: + Grid &grid; + Boundary &bc; + Velocities &velocities; + bool Steady{true}; + int iterations{-1}; + int innerIterations{1}; + T timestep{-1}; + int numThreads{omp_get_num_procs()}; + T porosity{1}; + const Eigen::MatrixX &outx; + const Eigen::MatrixX &outy; + CSV_OUTPUT csv_output{CSV_OUTPUT_OFF}; +}; +} // namespace tug diff --git a/include/tug/Core/Velocities.hpp b/include/tug/Core/Velocities.hpp new file mode 100644 index 0000000..51cc3a8 --- /dev/null +++ b/include/tug/Core/Velocities.hpp @@ -0,0 +1,488 @@ +/** + * @file Velocities.hpp + * @brief API of Velocities class, holding information for a simulation of + * Hydraulic Charge and Darcy-Velocities. Holds a predifined Grid object and + * Boundary object. + * + */ + +#pragma once + +#include "tug/Core/Matrix.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace Eigen; +namespace tug { +template class Velocities { +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) { + outx = MatrixX::Constant(grid.getRow(), grid.getCol() + 1, 0); + outy = MatrixX::Constant(grid.getRow() + 1, grid.getCol(), 0); + center = std::make_pair(grid.getRow() / 2, grid.getCol() / 2); + }; + + /** + * @brief Sets a fixed, constant hydraulic charge at domain center. + * + * @param inj_h fixed hydraulic charge at domain center. + */ + void setInjh(T inj_h) { + if (inj_h >= 0) { + this->inj_h = inj_h; + RowMajMat &concentrations = grid.getConcentrations(); + concentrations(center.first, center.second) = inj_h; + injhIsSet = true; + } else { + throw std::invalid_argument("Fixed hydraulic charge can not be negative"); + } + }; + + /** + * @brief Sets a constant permeability coefficient + * @param K constant permeability coefficient + */ + void setK(T K) { + this->K = K; + MatrixXd alphax = MatrixXd::Constant(grid.getRow(), grid.getCol(), K); + MatrixXd alphay = MatrixXd::Constant(grid.getRow(), grid.getCol(), K); + grid.setAlpha(alphax, alphay); + }; + + /** + * @brief Set the epsilon value, the relativ error allowed for convergence + * + * @param epsilon the new epsilon value + */ + void setEpsilon(T epsilon) { + if (0 <= epsilon && epsilon < 1) { + this->epsilon = epsilon; + } else { + throw std::invalid_argument( + "Relative Error epsilon must be between 0 and 1"); + } + } + + /** + * @brief Set the timestep per iteration + * + * @param timestep timestep per iteration + */ + void setTimestep(T timestep) { + 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 minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); + T cfl = minDeltaSquare / (4 * K); + if (timestep > cfl) { + this->innerIterations = (int)ceil(timestep / cfl); + this->timestep = timestep / (double)innerIterations; + std::cerr << "Warning :: Timestep was adjusted, because of stability " + "conditions. Time duration was approximately preserved by " + "adjusting internal number of iterations." + << std::endl; + std::cout << "FTCS" << "_" << "2D" << " :: Required " + << this->innerIterations + << " inner iterations with dt=" << this->timestep << std::endl; + } else { + this->innerIterations = 1; + } + }; + + /** + * @brief Set the number of iterations. If set to a number smaller than 1, + * calculation will terminate at convergence + * + * @param iterations Number of desired iterations + */ + + void setIterations(int iterations) { this->iterations = iterations; } + + /** + * @brief Set the number of desired openMP Threads. + * + * @param num_threads Number of desired threads. Must have a value between + * 1 and the maximum available number of processors. The + * maximum number of processors is set as the default case during Velocities + * construction. + */ + void setNumberThreads(int num_threads) { + if (num_threads > 0 && num_threads <= omp_get_num_procs()) { + this->numThreads = num_threads; + } else { + int maxThreadNumber = omp_get_num_procs(); + throw std::invalid_argument( + "Number of threads exceeds the number of processor cores (" + + std::to_string(maxThreadNumber) + ") or is less than 1."); + } + } + + /** + * @brief Set the option to output the results to a CSV file. Off by default. + * + * + * @param csv_output Valid output option. The following options can be set + * here: + * - CSV_OUTPUT_OFF: do not produce csv output + * - CSV_OUTPUT_ON: produce csv output with last + * charge and velocities matrizes + * - CSV_OUTPUT_VERBOSE: produce csv output with all + * charge matrizes and and last velocities matrix + */ + void setOutputCSV(CSV_OUTPUT csv_output) { + if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { + throw std::invalid_argument("Invalid CSV output option given!"); + } + this->csv_output = csv_output; + if (csv_output >= CSV_OUTPUT_ON) { + filename1 = createCSVfile("Charge"); + filename2 = createCSVfile("outx"); + filename3 = createCSVfile("outy"); + } + } + + /** + * @brief Creates a CSV file with a name containing the current simulation + * parameters. If the data name already exists, an additional counter + * is appended to the name. The name of the file is built up as follows: + * + + + + * +.csv + * + * @return string Filename with configured simulation parameters. + */ + std::string createCSVfile(std::string Type) const { + std::ofstream file; + int appendIdent = 0; + std::string appendIdentString; + + std::string row = std::to_string(grid.getRow()); + std::string col = std::to_string(grid.getCol()); + std::string numIterations = std::to_string(iterations); + + std::string filename = + Type + "_" + row + "_" + col + "_" + numIterations + ".csv"; + + while (std::filesystem::exists(filename)) { + appendIdent += 1; + appendIdentString = std::to_string(appendIdent); + filename = Type + "_" + row + "_" + col + "_" + numIterations + "-" + + appendIdentString + ".csv"; + } + + file.open(filename); + if (!file) { + exit(1); + } + + file.close(); + + return filename; + } + + /** + * @brief Writes the currently calculated Hydraulic Charge values of the grid + * into the CSV file with the passed filename. + * + * @param filename Name of the file to which the Hydraulic Charge values are + * to be written. + */ + void printChargeCSV(const std::string &filename) const { + std::ofstream file; + + file.open(filename, std::ios_base::app); + if (!file) { + exit(1); + } + + Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); + file << grid.getConcentrations().format(do_not_align) << std::endl; + file << std::endl << std::endl; + file.close(); + } + + /** + * @brief Reads a matrix stored in a CSV file and uses it for values of + * Hydraulic Heads, Matrix and grid must be of equal size + * + * @param filename name of the CSV file + */ + void readChargeCSV(std::string filename) { + std::ifstream file(filename); + std::string line; + Eigen::MatrixXd matrix; + + if (file.is_open()) { + while (std::getline(file, line)) { + std::istringstream iss(line); + double value; + std::vector row; + while (iss >> value) { + row.push_back(value); + } + if (!row.empty()) { + if (matrix.rows() == 0) { + matrix.resize(1, row.size()); + } else { + matrix.conservativeResize(matrix.rows() + 1, Eigen::NoChange); + } + matrix.row(matrix.rows() - 1) = + Eigen::VectorXd::Map(row.data(), row.size()); + } + } + file.close(); + } else { + std::cerr << "Unable to open file: " << filename << std::endl; + } + if (matrix.rows() == grid.getRow() && matrix.cols() == grid.getCol()) { + grid.setConcentrations(matrix); + velocities(); + if (csv_output > CSV_OUTPUT_OFF) { + printChargeCSV(filename1); + printVelocitiesCSV(filename2, filename3); + } + } else { + std::cerr << "gridsize and size of stored matrix dont align\n"; + } + } + + /** + * @brief Writes the current Darcy-velocities into a CSV file + * + * @param filenamex Name of the file to which velocities in direction x are + * written + * @param filenamey Name of the file to which velocities in direction y are + * written + */ + void printVelocitiesCSV(const std::string &filenamex, + const std::string &filenamey) const { + std::ofstream filex; + std::ofstream filey; + + filex.open(filenamex, std::ios_base::app); + if (!filex) { + exit(1); + } + filey.open(filenamey, std::ios_base::app); + if (!filey) { + exit(1); + } + + Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); + filex << outx.format(do_not_align) << std::endl; + filex << std::endl << std::endl; + filex.close(); + + filey << outy.format(do_not_align) << std::endl; + filey << std::endl << std::endl; + filey.close(); + } + + /** + * @brief Calculate the new hydraulic charge using FTCS + */ + void hydraulic_charge() { + FTCS_2D(this->grid, this->bc, this->timestep, this->numThreads); + if (injhIsSet == true) { + RowMajMat &concentrations = grid.getConcentrations(); + concentrations(center.first, center.second) = inj_h; + } + }; + + /** + * @brief checks if the matrix of Hydraulic Heads has converged + * + * @return bool true if for all corresponding cells of the matrices, + * containing old and new Charge values, the relative error is below the + * selected Epsilon + */ + bool checkConvergance(Eigen::MatrixX oldHeads, + Eigen::MatrixX newHeads) { + for (int i = 0; i < grid.getRow(); i++) { + for (int j = 0; j < grid.getCol(); j++) { + if (newHeads(i, j) != 0) { + if (abs((oldHeads(i, j) - newHeads(i, j)) / newHeads(i, j)) > + epsilon) { + return false; + } + } + } + } + return true; + } + + /** + * @brief Update the matrices containing Darcy velocities in x and y + * directions + */ + void velocities() { + int rows = grid.getRow(); + int cols = grid.getCol(); + float dx = grid.getDeltaRow(); + float dy = grid.getDeltaCol(); + Eigen::MatrixX concentrations = grid.getConcentrations(); +// calculate outx +#pragma omp parallel for num_threads(numThreads) + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols + 1; j++) { + if (j == 0) { + if (bc.getBoundaryElementType(BC_SIDE_LEFT, i) == BC_TYPE_CLOSED) { + outx(i, j) = 0; + } else { + outx(i, j) = -K * + (concentrations(i, j) - + bc.getBoundaryElementValue(BC_SIDE_LEFT, i)) / + (dx / 2); + } + } else if (j == cols) { + if (bc.getBoundaryElementType(BC_SIDE_RIGHT, i) == BC_TYPE_CLOSED) { + outx(i, j) = 0; + } else { + outx(i, j) = -K * + (bc.getBoundaryElementValue(BC_SIDE_RIGHT, i) - + concentrations(i, j - 1)) / + (dx / 2); + } + + } else { + outx(i, j) = + -K * (concentrations(i, j) - concentrations(i, j - 1)) / dx; + } + } + } +// calculate outy +#pragma omp parallel for num_threads(numThreads) + for (int i = 0; i < rows + 1; i++) { + for (int j = 0; j < cols; j++) { + if (i == 0) { + if (bc.getBoundaryElementType(BC_SIDE_TOP, j) == BC_TYPE_CLOSED) { + outy(i, j) = 0; + } else { + outy(i, j) = -K * + (concentrations(i, j) - + bc.getBoundaryElementValue(BC_SIDE_TOP, j)) / + (dy / 2); + } + } else if (i == rows) { + if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, j) == BC_TYPE_CLOSED) { + outy(i, j) = 0; + } else { + outy(i, j) = -K * + (bc.getBoundaryElementValue(BC_SIDE_BOTTOM, j) - + concentrations(i - 1, j)) / + (dy / 2); + } + } else { + outy(i, j) = + -K * (concentrations(i, j) - concentrations(i - 1, j)) / dy; + } + } + } + }; + + /** + * @brief Getter function for outx, the matrix containing velocities in + * x-Direction; returns a reference to outx + * + * */ + const Eigen::MatrixX &getOutx() const { return outx; } + + /** + * @brief Getter function for outy, the matrix containing velocities in + * y-Direction; return a reference to outy + */ + const Eigen::MatrixX &getOuty() const { return outy; } + + /** + * @brief Simulation of hydraulic charge either until convergence, + * or for a number of selected timesteps. Calculation of Darcy-velocities. + */ + void run() { + // if iterations < 1 calculate hydraulic charge until steady state is + // reached + if (iterations < 1) { + // Calculate largest possible timestep, depending on K and gridsize + const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); + const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); + const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); + setTimestep(minDeltaSquare / (4 * K)); + + Eigen::MatrixX oldConcentrations; + do { + oldConcentrations = grid.getConcentrations(); + if (csv_output >= CSV_OUTPUT_VERBOSE) { + printChargeCSV(filename1); + } + for (int i = 0; i < (grid.getRow() + grid.getCol() - 2); i++) { + hydraulic_charge(); + } + } while (checkConvergance(oldConcentrations, grid.getConcentrations()) == + false); + } + // if iterations >= 1 calculate hydraulice charge for a given number of + // iterations + else { + if (timestep == -1) { + throw_invalid_argument("Timestep is not set"); + } + for (int i = 0; i < iterations * innerIterations; i++) { + hydraulic_charge(); + } + } + + velocities(); + + if (csv_output > CSV_OUTPUT_OFF) { + printChargeCSV(filename1); + printVelocitiesCSV(filename2, filename3); + } + }; + +private: + int iterations{-1}; + int innerIterations{1}; + bool injhIsSet{false}; + T timestep{-1}; + T inj_h{1}; + T K{1}; + T epsilon{1E-5}; + int numThreads{omp_get_num_procs()}; + Grid &grid; + Boundary &bc; + CSV_OUTPUT csv_output{CSV_OUTPUT_OFF}; + Eigen::MatrixX outx; + Eigen::MatrixX outy; + std::pair center; + std::string filename1; + std::string filename2; + std::string filename3; +}; +} // namespace tug diff --git a/include/tug/Matrix.hpp b/include/tug/Matrix.hpp new file mode 100644 index 0000000..e69de29 From 3b953e0b96caa0899fc38cf9e8d12f19e664277f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 10 Dec 2024 12:00:24 +0100 Subject: [PATCH 02/20] refactor: Consolidate includes by introducing tug.hpp for cleaner code structure --- examples/Advection.cpp | 3 +-- examples/BTCS_2D_proto_example.cpp | 2 +- examples/FTCS_2D_proto_closed_mdl.cpp | 2 +- examples/FTCS_2D_proto_example_mdl.cpp | 2 +- include/tug/{ => Advection}/Advection.hpp | 4 ++-- include/tug/{Core => Advection}/Velocities.hpp | 2 +- include/tug/{ => Diffusion}/Diffusion.hpp | 0 include/tug/Matrix.hpp | 0 include/tug/tug.hpp | 7 +++++++ test/testDiffusion.cpp | 3 ++- 10 files changed, 16 insertions(+), 9 deletions(-) rename include/tug/{ => Advection}/Advection.hpp (99%) rename include/tug/{Core => Advection}/Velocities.hpp (99%) rename include/tug/{ => Diffusion}/Diffusion.hpp (100%) delete mode 100644 include/tug/Matrix.hpp create mode 100644 include/tug/tug.hpp diff --git a/examples/Advection.cpp b/examples/Advection.cpp index 414a8af..7876243 100644 --- a/examples/Advection.cpp +++ b/examples/Advection.cpp @@ -1,7 +1,6 @@ #include #include -#include -#include +#include using namespace Eigen; using namespace tug; diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index c09d93d..f793e5c 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -1,5 +1,5 @@ #include -#include +#include using namespace Eigen; using namespace tug; diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index 69b31c3..eab6dc9 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -9,7 +9,7 @@ #include #include #include -#include +#include using namespace Eigen; using namespace tug; diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index 41ddb91..f7059a1 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -7,7 +7,7 @@ */ #include -#include +#include using namespace Eigen; using namespace tug; diff --git a/include/tug/Advection.hpp b/include/tug/Advection/Advection.hpp similarity index 99% rename from include/tug/Advection.hpp rename to include/tug/Advection/Advection.hpp index 4d530bf..855e1b5 100644 --- a/include/tug/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -23,9 +23,9 @@ #include #include #include -#include +#include -#include +#include using namespace Eigen; namespace tug { diff --git a/include/tug/Core/Velocities.hpp b/include/tug/Advection/Velocities.hpp similarity index 99% rename from include/tug/Core/Velocities.hpp rename to include/tug/Advection/Velocities.hpp index 51cc3a8..dbfef09 100644 --- a/include/tug/Core/Velocities.hpp +++ b/include/tug/Advection/Velocities.hpp @@ -25,7 +25,7 @@ #include #include #include -#include +#include using namespace Eigen; namespace tug { diff --git a/include/tug/Diffusion.hpp b/include/tug/Diffusion/Diffusion.hpp similarity index 100% rename from include/tug/Diffusion.hpp rename to include/tug/Diffusion/Diffusion.hpp diff --git a/include/tug/Matrix.hpp b/include/tug/Matrix.hpp deleted file mode 100644 index e69de29..0000000 diff --git a/include/tug/tug.hpp b/include/tug/tug.hpp new file mode 100644 index 0000000..a766c3a --- /dev/null +++ b/include/tug/tug.hpp @@ -0,0 +1,7 @@ +#pragma once + +#include +#include +#include +#include +#include \ No newline at end of file diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 15fa18b..c49a166 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -3,7 +3,8 @@ #include "gtest/gtest.h" #include #include -#include + +#include #include #include From 031c1b2eef27668b130b48d1ad6e24498ff4e62d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 11 Dec 2024 13:22:23 +0100 Subject: [PATCH 03/20] feat: Implement inner boundaries for FTCS --- include/tug/Core/Numeric/FTCS.hpp | 24 ++++++++++++++++++++++++ test/testDiffusion.cpp | 30 ++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+) diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp index ce1d5e8..3ab5cc2 100644 --- a/include/tug/Core/Numeric/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -11,6 +11,7 @@ #include "tug/Core/TugUtils.hpp" #include #include +#include #include #include #include @@ -55,6 +56,27 @@ constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center, tug_assert(false, "Undefined Boundary Condition Type!"); } +template +static inline void checkAndSetConstantInnerCells(const Boundary &bc, + Grid &grid) { + 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; + concentrations(row, col) = value; + } + + return; +} + // FTCS solution for 1D grid template static void FTCS_1D(SimulationInput &input) { const std::size_t &colMax = input.colMax; @@ -71,6 +93,8 @@ template static void FTCS_1D(SimulationInput &input) { // only one row in 1D case -> row constant at index 0 int row = 0; + const auto inner_bc = bc.getInnerBoundaryRow(0); + // inner cells // independent of boundary condition type for (int col = 1; col < colMax - 1; col++) { diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index c49a166..9a74670 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -215,3 +215,33 @@ 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 From 13d55f9260cc9d8c5c25dbde561ac2a664db3d6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 11 Dec 2024 15:30:02 +0100 Subject: [PATCH 04/20] refactor: Velocities.hpp --- include/tug/Advection/Velocities.hpp | 519 +++++++++------------------ 1 file changed, 169 insertions(+), 350 deletions(-) diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp index dbfef09..dd495fa 100644 --- a/include/tug/Advection/Velocities.hpp +++ b/include/tug/Advection/Velocities.hpp @@ -8,28 +8,47 @@ #pragma once -#include "tug/Core/Matrix.hpp" #include -#include -#include +#include #include #include #include #include -#include -#include -#include -#include -#include +#include +#include #include #include #include #include +#include + +#ifdef _OPENMP +#include +#else +#define omp_get_num_procs() 1 +#endif using namespace Eigen; namespace tug { -template class Velocities { + +enum HYDRAULIC_MODE { TRANSIENT, STEADY_STATE }; +enum HYDRAULIC_RESOLVE { EXPLICIT, IMPLICIT }; + +template +class Velocities : public BaseSimulation { +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; + public: /** * @brief Construct a new Velocities object, used to calculate Hydraulic @@ -42,38 +61,9 @@ public: * @param grid Valid grid object * @param bc Valid boundary condition object */ - Velocities(Grid &_grid, Boundary &_bc) : grid(_grid), bc(_bc) { - outx = MatrixX::Constant(grid.getRow(), grid.getCol() + 1, 0); - outy = MatrixX::Constant(grid.getRow() + 1, grid.getCol(), 0); - center = std::make_pair(grid.getRow() / 2, grid.getCol() / 2); - }; - - /** - * @brief Sets a fixed, constant hydraulic charge at domain center. - * - * @param inj_h fixed hydraulic charge at domain center. - */ - void setInjh(T inj_h) { - if (inj_h >= 0) { - this->inj_h = inj_h; - RowMajMat &concentrations = grid.getConcentrations(); - concentrations(center.first, center.second) = inj_h; - injhIsSet = true; - } else { - throw std::invalid_argument("Fixed hydraulic charge can not be negative"); - } - }; - - /** - * @brief Sets a constant permeability coefficient - * @param K constant permeability coefficient - */ - void setK(T K) { - this->K = K; - MatrixXd alphax = MatrixXd::Constant(grid.getRow(), grid.getCol(), K); - MatrixXd alphay = MatrixXd::Constant(grid.getRow(), grid.getCol(), K); - grid.setAlpha(alphax, alphay); - }; + Velocities(Grid &_grid, Boundary &_bc) + : grid(_grid), bc(_bc), velocitiesX(grid.getRow(), grid.getCol() + 1), + velocitiesY(grid.getRow() + 1, grid.getCol()) {}; /** * @brief Set the epsilon value, the relativ error allowed for convergence @@ -81,12 +71,10 @@ public: * @param epsilon the new epsilon value */ void setEpsilon(T epsilon) { - if (0 <= epsilon && epsilon < 1) { - this->epsilon = epsilon; - } else { - throw std::invalid_argument( - "Relative Error epsilon must be between 0 and 1"); - } + tug_assert(0 <= epsilon && epsilon < 1, + "Relative Error epsilon must be between 0 and 1"); + + this->epsilon = epsilon; } /** @@ -102,7 +90,11 @@ public: const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); - T cfl = minDeltaSquare / (4 * K); + + const T maxK = + std::max(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff()); + + T cfl = minDeltaSquare / (4 * maxK); if (timestep > cfl) { this->innerIterations = (int)ceil(timestep / cfl); this->timestep = timestep / (double)innerIterations; @@ -118,15 +110,6 @@ public: } }; - /** - * @brief Set the number of iterations. If set to a number smaller than 1, - * calculation will terminate at convergence - * - * @param iterations Number of desired iterations - */ - - void setIterations(int iterations) { this->iterations = iterations; } - /** * @brief Set the number of desired openMP Threads. * @@ -147,173 +130,63 @@ public: } /** - * @brief Set the option to output the results to a CSV file. Off by default. + * @brief Getter function for outx, the matrix containing velocities in + * x-Direction; returns a reference to outx * - * - * @param csv_output Valid output option. The following options can be set - * here: - * - CSV_OUTPUT_OFF: do not produce csv output - * - CSV_OUTPUT_ON: produce csv output with last - * charge and velocities matrizes - * - CSV_OUTPUT_VERBOSE: produce csv output with all - * charge matrizes and and last velocities matrix - */ - void setOutputCSV(CSV_OUTPUT csv_output) { - if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { - throw std::invalid_argument("Invalid CSV output option given!"); - } - this->csv_output = csv_output; - if (csv_output >= CSV_OUTPUT_ON) { - filename1 = createCSVfile("Charge"); - filename2 = createCSVfile("outx"); - filename3 = createCSVfile("outy"); - } - } + * */ + const RowMajMat &getVelocitiesX() const { return this->velocitiesX; } /** - * @brief Creates a CSV file with a name containing the current simulation - * parameters. If the data name already exists, an additional counter - * is appended to the name. The name of the file is built up as follows: - * + + + - * +.csv - * - * @return string Filename with configured simulation parameters. + * @brief Getter function for outy, the matrix containing velocities in + * y-Direction; return a reference to outy */ - std::string createCSVfile(std::string Type) const { - std::ofstream file; - int appendIdent = 0; - std::string appendIdentString; - - std::string row = std::to_string(grid.getRow()); - std::string col = std::to_string(grid.getCol()); - std::string numIterations = std::to_string(iterations); - - std::string filename = - Type + "_" + row + "_" + col + "_" + numIterations + ".csv"; - - while (std::filesystem::exists(filename)) { - appendIdent += 1; - appendIdentString = std::to_string(appendIdent); - filename = Type + "_" + row + "_" + col + "_" + numIterations + "-" + - appendIdentString + ".csv"; - } - - file.open(filename); - if (!file) { - exit(1); - } - - file.close(); - - return filename; - } + const RowMajMat &getVelocitiesY() const { return this->velocitiesY; } /** - * @brief Writes the currently calculated Hydraulic Charge values of the grid - * into the CSV file with the passed filename. - * - * @param filename Name of the file to which the Hydraulic Charge values are - * to be written. + * @brief Simulation of hydraulic charge either until convergence, + * or for a number of selected timesteps. Calculation of Darcy-velocities. */ - void printChargeCSV(const std::string &filename) const { - std::ofstream file; + void run() { + // if iterations < 1 calculate hydraulic charge until steady state is + // reached - file.open(filename, std::ios_base::app); - if (!file) { - exit(1); - } + 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(); + const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); - Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); - file << grid.getConcentrations().format(do_not_align) << std::endl; - file << std::endl << std::endl; - file.close(); - } + const T maxK = + std::max(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff()); - /** - * @brief Reads a matrix stored in a CSV file and uses it for values of - * Hydraulic Heads, Matrix and grid must be of equal size - * - * @param filename name of the CSV file - */ - void readChargeCSV(std::string filename) { - std::ifstream file(filename); - std::string line; - Eigen::MatrixXd matrix; + setTimestep(minDeltaSquare / (4 * maxK)); - if (file.is_open()) { - while (std::getline(file, line)) { - std::istringstream iss(line); - double value; - std::vector row; - while (iss >> value) { - row.push_back(value); - } - if (!row.empty()) { - if (matrix.rows() == 0) { - matrix.resize(1, row.size()); - } else { - matrix.conservativeResize(matrix.rows() + 1, Eigen::NoChange); - } - matrix.row(matrix.rows() - 1) = - Eigen::VectorXd::Map(row.data(), row.size()); - } - } - file.close(); + RowMajMat oldConcentrations; + do { + oldConcentrations = grid.getConcentrations(); + (void)calculate_hydraulic_flow(); + } while (!checkConvergance(oldConcentrations, grid.getConcentrations())); } else { - std::cerr << "Unable to open file: " << filename << std::endl; - } - if (matrix.rows() == grid.getRow() && matrix.cols() == grid.getCol()) { - grid.setConcentrations(matrix); - velocities(); - if (csv_output > CSV_OUTPUT_OFF) { - printChargeCSV(filename1); - printVelocitiesCSV(filename2, filename3); + if (timestep == -1) { + throw_invalid_argument("Timestep is not set"); + } + for (int i = 0; i < innerIterations; i++) { + (void)calculate_hydraulic_flow(); } - } else { - std::cerr << "gridsize and size of stored matrix dont align\n"; - } - } - - /** - * @brief Writes the current Darcy-velocities into a CSV file - * - * @param filenamex Name of the file to which velocities in direction x are - * written - * @param filenamey Name of the file to which velocities in direction y are - * written - */ - void printVelocitiesCSV(const std::string &filenamex, - const std::string &filenamey) const { - std::ofstream filex; - std::ofstream filey; - - filex.open(filenamex, std::ios_base::app); - if (!filex) { - exit(1); - } - filey.open(filenamey, std::ios_base::app); - if (!filey) { - exit(1); } - Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); - filex << outx.format(do_not_align) << std::endl; - filex << std::endl << std::endl; - filex.close(); - - filey << outy.format(do_not_align) << std::endl; - filey << std::endl << std::endl; - filey.close(); - } + (void)computeFluidVelocities(); + }; +private: /** * @brief Calculate the new hydraulic charge using FTCS */ - void hydraulic_charge() { - FTCS_2D(this->grid, this->bc, this->timestep, this->numThreads); - if (injhIsSet == true) { - RowMajMat &concentrations = grid.getConcentrations(); - concentrations(center.first, center.second) = inj_h; + void calculate_hydraulic_flow() { + if constexpr (hyd_resolve == EXPLICIT) { + FTCS_2D(this->grid, this->bc, this->timestep, this->numThreads); + } else { + BTCS_2D(this->grid, this->bc, this->timestep, ThomasAlgorithm); } }; @@ -324,165 +197,111 @@ public: * containing old and new Charge values, the relative error is below the * selected Epsilon */ - bool checkConvergance(Eigen::MatrixX oldHeads, - Eigen::MatrixX newHeads) { - for (int i = 0; i < grid.getRow(); i++) { - for (int j = 0; j < grid.getCol(); j++) { - if (newHeads(i, j) != 0) { - if (abs((oldHeads(i, j) - newHeads(i, j)) / newHeads(i, j)) > - epsilon) { - return false; - } - } - } - } - return true; + bool checkConvergance(const RowMajMat &oldHeads, + const RowMajMat &newHeads) { + const auto abs_err = (oldHeads - newHeads).cwiseAbs(); + const auto rel_err = abs_err.cwiseQuotient(newHeads); + + return rel_err.maxCoeff() < epsilon; } /** * @brief Update the matrices containing Darcy velocities in x and y * directions */ - void velocities() { - int rows = grid.getRow(); - int cols = grid.getCol(); - float dx = grid.getDeltaRow(); - float dy = grid.getDeltaCol(); - Eigen::MatrixX concentrations = grid.getConcentrations(); -// calculate outx + 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 RowMajMat &permKX = grid.getAlphaX(); + const RowMajMat &permKY = grid.getAlphaY(); + + // calculate velocities in x-direction + for (std::size_t i_rows = 0; i_rows < rows; i_rows++) { + const auto bc_left = bc.getBoundaryElement(BC_SIDE_LEFT, i_rows); + switch (bc_left.getType()) { + case BC_TYPE_CLOSED: { + velocitiesX(i_rows, 0) = 0; + break; + } + case BC_TYPE_CONSTANT: { + velocitiesX(i_rows, 0) = + -permKX(i_rows, 0) * + (hydraulicCharges(i_rows, 0) - bc_left.getValue()) / (dx / 2); + break; + } + } + + const auto bc_right = bc.getBoundaryElement(BC_SIDE_RIGHT, i_rows); + switch (bc_right.getType()) { + case BC_TYPE_CLOSED: { + velocitiesX(i_rows, cols) = 0; + break; + } + case BC_TYPE_CONSTANT: { + velocitiesX(i_rows, cols) = + -permKX(i_rows, cols - 1) * + (bc_right.getValue() - hydraulicCharges(i_rows, cols - 1)) / + (dx / 2); + break; + } + } + } + +// main loop for calculating velocities in x-direction for inner cells #pragma omp parallel for num_threads(numThreads) for (int i = 0; i < rows; i++) { - for (int j = 0; j < cols + 1; j++) { - if (j == 0) { - if (bc.getBoundaryElementType(BC_SIDE_LEFT, i) == BC_TYPE_CLOSED) { - outx(i, j) = 0; - } else { - outx(i, j) = -K * - (concentrations(i, j) - - bc.getBoundaryElementValue(BC_SIDE_LEFT, i)) / - (dx / 2); - } - } else if (j == cols) { - if (bc.getBoundaryElementType(BC_SIDE_RIGHT, i) == BC_TYPE_CLOSED) { - outx(i, j) = 0; - } else { - outx(i, j) = -K * - (bc.getBoundaryElementValue(BC_SIDE_RIGHT, i) - - concentrations(i, j - 1)) / - (dx / 2); - } - - } else { - outx(i, j) = - -K * (concentrations(i, j) - concentrations(i, j - 1)) / dx; - } + for (int j = 1; j < cols; j++) { + velocitiesX(i, j) = + -permKX(i, j - 1) * + (hydraulicCharges(i, j) - hydraulicCharges(i, j - 1)) / dx; } } -// calculate outy + + // calculate velocities in y-direction + for (std::size_t i_cols = 0; i_cols < cols; i_cols++) { + const auto bc_top = bc.getBoundaryElement(BC_SIDE_TOP, i_cols); + switch (bc_top.getType()) { + case BC_TYPE_CLOSED: { + velocitiesY(0, i_cols) = 0; + break; + } + case BC_TYPE_CONSTANT: { + velocitiesY(0, i_cols) = + -permKY(0, i_cols) * + (hydraulicCharges(0, i_cols) - bc_top.getValue()) / (dy / 2); + break; + } + } + + const auto bc_bottom = bc.getBoundaryElement(BC_SIDE_BOTTOM, i_cols); + switch (bc_bottom.getType()) { + case BC_TYPE_CLOSED: { + velocitiesY(rows, i_cols) = 0; + break; + } + case BC_TYPE_CONSTANT: { + velocitiesY(rows, i_cols) = + -permKY(rows - 1, i_cols) * + (bc_bottom.getValue() - hydraulicCharges(rows - 1, i_cols)) / + (dy / 2); + break; + } + } + } + +// main loop for calculating velocities in y-direction for inner cells #pragma omp parallel for num_threads(numThreads) - for (int i = 0; i < rows + 1; i++) { + for (int i = 1; i < rows; i++) { for (int j = 0; j < cols; j++) { - if (i == 0) { - if (bc.getBoundaryElementType(BC_SIDE_TOP, j) == BC_TYPE_CLOSED) { - outy(i, j) = 0; - } else { - outy(i, j) = -K * - (concentrations(i, j) - - bc.getBoundaryElementValue(BC_SIDE_TOP, j)) / - (dy / 2); - } - } else if (i == rows) { - if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, j) == BC_TYPE_CLOSED) { - outy(i, j) = 0; - } else { - outy(i, j) = -K * - (bc.getBoundaryElementValue(BC_SIDE_BOTTOM, j) - - concentrations(i - 1, j)) / - (dy / 2); - } - } else { - outy(i, j) = - -K * (concentrations(i, j) - concentrations(i - 1, j)) / dy; - } + velocitiesY(i, j) = + -permKY(i, j) * + (hydraulicCharges(i, j) - hydraulicCharges(i - 1, j)) / dy; } } }; - - /** - * @brief Getter function for outx, the matrix containing velocities in - * x-Direction; returns a reference to outx - * - * */ - const Eigen::MatrixX &getOutx() const { return outx; } - - /** - * @brief Getter function for outy, the matrix containing velocities in - * y-Direction; return a reference to outy - */ - const Eigen::MatrixX &getOuty() const { return outy; } - - /** - * @brief Simulation of hydraulic charge either until convergence, - * or for a number of selected timesteps. Calculation of Darcy-velocities. - */ - void run() { - // if iterations < 1 calculate hydraulic charge until steady state is - // reached - if (iterations < 1) { - // Calculate largest possible timestep, depending on K and gridsize - const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); - const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); - const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); - setTimestep(minDeltaSquare / (4 * K)); - - Eigen::MatrixX oldConcentrations; - do { - oldConcentrations = grid.getConcentrations(); - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printChargeCSV(filename1); - } - for (int i = 0; i < (grid.getRow() + grid.getCol() - 2); i++) { - hydraulic_charge(); - } - } while (checkConvergance(oldConcentrations, grid.getConcentrations()) == - false); - } - // if iterations >= 1 calculate hydraulice charge for a given number of - // iterations - else { - if (timestep == -1) { - throw_invalid_argument("Timestep is not set"); - } - for (int i = 0; i < iterations * innerIterations; i++) { - hydraulic_charge(); - } - } - - velocities(); - - if (csv_output > CSV_OUTPUT_OFF) { - printChargeCSV(filename1); - printVelocitiesCSV(filename2, filename3); - } - }; - -private: - int iterations{-1}; - int innerIterations{1}; - bool injhIsSet{false}; - T timestep{-1}; - T inj_h{1}; - T K{1}; - T epsilon{1E-5}; - int numThreads{omp_get_num_procs()}; - Grid &grid; - Boundary &bc; - CSV_OUTPUT csv_output{CSV_OUTPUT_OFF}; - Eigen::MatrixX outx; - Eigen::MatrixX outy; - std::pair center; - std::string filename1; - std::string filename2; - std::string filename3; }; } // namespace tug From 1a11991af0b52a58011569f3021753c3fa075696 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 11 Dec 2024 15:30:10 +0100 Subject: [PATCH 05/20] feat: Add unit tests for Velocities functionality --- include/tug/tug.hpp | 2 +- test/CMakeLists.txt | 1 + test/testVelocities.cpp | 84 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 86 insertions(+), 1 deletion(-) create mode 100644 test/testVelocities.cpp diff --git a/include/tug/tug.hpp b/include/tug/tug.hpp index a766c3a..89440b9 100644 --- a/include/tug/tug.hpp +++ b/include/tug/tug.hpp @@ -1,6 +1,6 @@ #pragma once -#include +// #include #include #include #include diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 2a8d780..a932356 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -14,6 +14,7 @@ setup.cpp testDiffusion.cpp testFTCS.cpp testBoundary.cpp +testVelocities.cpp ) target_link_libraries(testTug tug GTest::gtest) diff --git a/test/testVelocities.cpp b/test/testVelocities.cpp new file mode 100644 index 0000000..03b1450 --- /dev/null +++ b/test/testVelocities.cpp @@ -0,0 +1,84 @@ +#include "tug/Boundary.hpp" +#include "tug/Core/Matrix.hpp" +#include +#include + +#include + +#define VELOCITIES_TEST(x) TEST(Velocities, x) + +VELOCITIES_TEST(SteadyStateCenter) { + // Arrange + constexpr std::size_t rows = 25; + constexpr std::size_t cols = 25; + + constexpr std::size_t center_row = rows / 2; + constexpr std::size_t center_col = cols / 2; + + constexpr double K = 1E-2; + + tug::RowMajMat hydHeads = + tug::RowMajMat::Constant(rows, cols, 1); + + tug::RowMajMat permKX = + tug::RowMajMat::Constant(rows, cols, K); + tug::RowMajMat permKY = + tug::RowMajMat::Constant(rows, cols, K); + + tug::Grid64 gridHeads(hydHeads); + gridHeads.setDomain(100, 100); + gridHeads.setAlpha(permKX, permKY); + + tug::Boundary bcH = tug::Boundary(gridHeads); + bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1); + bcH.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1); + bcH.setBoundarySideConstant(tug::BC_SIDE_TOP, 1); + bcH.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1); + + bcH.setInnerBoundary(center_row, center_col, 10); + + tug::Velocities + velocities(gridHeads, bcH); + + velocities.run(); + + const auto &velocitiesX = velocities.getVelocitiesX(); + const auto &velocitiesY = velocities.getVelocitiesY(); + + // Assert + + // check velocities in x-direction + for (std::size_t i_rows = 0; i_rows < rows; i_rows++) { + for (std::size_t i_cols = 0; i_cols < cols + 1; i_cols++) { + if (i_rows <= center_row && i_cols <= center_col) { + EXPECT_LE(velocitiesX(i_rows, i_cols), 0); + } else if (i_rows > center_row && i_cols > center_col) { + EXPECT_GE(velocitiesX(i_rows, i_cols), 0); + } else if (i_rows <= center_row && i_cols > center_col) { + EXPECT_GE(velocitiesX(i_rows, i_cols), 0); + } else if (i_rows > center_row && i_cols <= center_col) { + EXPECT_LE(velocitiesX(i_rows, i_cols), 0); + } else { + FAIL() << "Uncovered case"; + } + } + } + + // check velocities in y-direction + for (std::size_t i_rows = 0; i_rows < rows + 1; i_rows++) { + for (std::size_t i_cols = 0; i_cols < cols; i_cols++) { + if (i_rows <= center_row && i_cols <= center_col) { + EXPECT_LE(velocitiesY(i_rows, i_cols), 0); + } else if (i_rows > center_row && i_cols > center_col) { + EXPECT_GE(velocitiesY(i_rows, i_cols), 0); + } else if (i_rows <= center_row && i_cols > center_col) { + EXPECT_LE(velocitiesY(i_rows, i_cols), 0); + } else if (i_rows > center_row && i_cols <= center_col) { + EXPECT_GE(velocitiesY(i_rows, i_cols), 0); + } else { + FAIL() << "Uncovered case"; + } + } + } +} \ No newline at end of file From ca94cebba290c350d7a160aee7766259e5651765 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 5 Feb 2025 13:53:08 +0100 Subject: [PATCH 06/20] chore: add missing cstdint include --- include/tug/Core/Numeric/SimulationInput.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/tug/Core/Numeric/SimulationInput.hpp b/include/tug/Core/Numeric/SimulationInput.hpp index 202619c..d68effc 100644 --- a/include/tug/Core/Numeric/SimulationInput.hpp +++ b/include/tug/Core/Numeric/SimulationInput.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -18,4 +19,4 @@ template struct SimulationInput { const T deltaRow; const T deltaCol; }; -} // namespace tug \ No newline at end of file +} // namespace tug 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 07/20] 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 +} From b7fcfb3ca52f4b673b34afb1c1152f5ec33920b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 6 Feb 2025 12:04:32 +0100 Subject: [PATCH 08/20] refactor(advection): rename alpha to permK for permeability --- include/tug/Advection/Velocities.hpp | 42 ++++++++++++++-------------- test/testVelocities.cpp | 4 +-- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp index 2fbd774..8cf2fc3 100644 --- a/include/tug/Advection/Velocities.hpp +++ b/include/tug/Advection/Velocities.hpp @@ -47,24 +47,24 @@ private: RowMajMat velocitiesX; RowMajMat velocitiesY; - RowMajMat alphaX; - RowMajMat alphaY; + RowMajMat permKX; + RowMajMat permKY; public: 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()) {}; + permKX(origin.rows(), origin.cols()), + permKY(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) {}; + velocitiesY(rows + 1, cols), permKX(rows, cols), permKY(rows, cols) {}; - Velocities(T *data, std::size_t length) - : BaseSimulationGrid(data, 1, length), velocitiesX(1, length + 1), - alphaX(1, length) {}; + // 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 @@ -83,19 +83,19 @@ public: * * @return RowMajMat& Reference to the alphaX matrix. */ - RowMajMat &getAlphaX() { return alphaX; } + RowMajMat &getPermKX() { return permKX; } /** * @brief Get the alphaY matrix. * * @return RowMajMat& Reference to the alphaY matrix. */ - RowMajMat &getAlphaY() { + RowMajMat &getPermKY() { tug_assert( this->getDim(), "Grid is not two dimensional, there is no domain in y-direction!"); - return alphaY; + return permKY; } /** @@ -103,19 +103,19 @@ public: * * @param alphaX The new alphaX matrix. */ - void setAlphaX(const RowMajMat &alphaX) { this->alphaX = alphaX; } + void setPermKX(const RowMajMat &alphaX) { this->permKX = alphaX; } /** * @brief Set the alphaY matrix. * * @param alphaY The new alphaY matrix. */ - void setAlphaY(const RowMajMat &alphaY) { + void setPermKY(const RowMajMat &alphaY) { tug_assert( this->getDim(), "Grid is not two dimensional, there is no domain in y-direction!"); - this->alphaY = alphaY; + this->permKY = alphaY; } /** @@ -132,7 +132,7 @@ public: const T deltaRowSquare = this->deltaRow() * this->deltaRow(); const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); - const T maxK = std::max(this->alphaX.maxCoeff(), this->alphaY.maxCoeff()); + const T maxK = std::max(this->permKX.maxCoeff(), this->permKY.maxCoeff()); T cfl = minDeltaSquare / (4 * maxK); if (timestep > cfl) { @@ -192,8 +192,8 @@ public: SimulationInput input = {.concentrations = this->getConcentrationMatrix(), - .alphaX = this->getAlphaX(), - .alphaY = this->getAlphaY(), + .alphaX = this->getPermKX(), + .alphaY = this->getPermKY(), .boundaries = this->getBoundaryConditions(), .dim = this->getDim(), .timestep = this->timestep, @@ -207,7 +207,7 @@ public: const T deltaRowSquare = this->deltaRow() * this->deltaRow(); const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); - const T maxK = std::max(this->alphaX.maxCoeff(), this->alphaY.maxCoeff()); + const T maxK = std::max(this->permKX.maxCoeff(), this->permKY.maxCoeff()); // Calculate largest possible timestep, depending on K and gridsize setTimestep(minDeltaSquare / (4 * maxK)); @@ -268,8 +268,8 @@ private: const T dy = this->deltaRow(); const RowMajMat &hydraulicCharges = this->getConcentrationMatrix(); - const RowMajMat &permKX = this->alphaX; - const RowMajMat &permKY = this->alphaY; + const RowMajMat &permKX = this->permKX; + const RowMajMat &permKY = this->permKY; const Boundary &bc = this->getBoundaryConditions(); @@ -283,7 +283,7 @@ private: } case BC_TYPE_CONSTANT: { velocitiesX(i_rows, 0) = - -this->alphaX(i_rows, 0) * + -this->permKX(i_rows, 0) * (hydraulicCharges(i_rows, 0) - bc_left.getValue()) / (dx / 2); break; } diff --git a/test/testVelocities.cpp b/test/testVelocities.cpp index 9e4ecac..0c8af66 100644 --- a/test/testVelocities.cpp +++ b/test/testVelocities.cpp @@ -30,8 +30,8 @@ VELOCITIES_TEST(SteadyStateCenter) { velo(hydHeads); velo.setDomain(100, 100); - velo.setAlphaX(permKX); - velo.setAlphaY(permKY); + velo.setPermKX(permKX); + velo.setPermKY(permKY); tug::Boundary &bcH = velo.getBoundaryConditions(); bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1); From 1ca81b440600e4c0f38939d3710bffae38949cbe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 6 Feb 2025 13:00:28 +0100 Subject: [PATCH 09/20] feat: Implement advection simulation with velocities and boundary conditions There is a bug that gains concentration even when inflow=outflow --- examples/Advection.cpp | 69 ++++---- examples/CMakeLists.txt | 12 +- include/tug/Advection/Advection.hpp | 249 +++++++++------------------- include/tug/Core/Numeric/FTCS.hpp | 31 ++-- include/tug/tug.hpp | 3 +- 5 files changed, 137 insertions(+), 227 deletions(-) diff --git a/examples/Advection.cpp b/examples/Advection.cpp index 7876243..c2ca66a 100644 --- a/examples/Advection.cpp +++ b/examples/Advection.cpp @@ -1,3 +1,6 @@ +#include "tug/Advection/Advection.hpp" +#include "tug/Advection/Velocities.hpp" +#include "tug/Core/Matrix.hpp" #include #include #include @@ -6,53 +9,53 @@ using namespace Eigen; using namespace tug; int main(int argc, char *argv[]) { - int row = 21; - int col = 21; + int row = 5; + int col = 5; // create two grids of equal size, grid1 for hydraulics heads, grid2 for // Concentrations RowMajMat hydHeads = RowMajMat::Constant(row, col, 1); RowMajMat concentrations = RowMajMat::Constant(row, col, 0); + hydHeads(row / 2, col / 2) = 10; + concentrations(row / 2, col / 2) = 1; - Grid64 gridHeads(hydHeads); - Grid64 gridConc(concentrations); - gridHeads.setDomain(100, 100); - gridConc.setDomain(100, 100); + Velocities + velocities(hydHeads); + velocities.setDomain(100, 100); + velocities.setPermKX(RowMajMat::Constant(row, col, 1E-8)); + velocities.setPermKY(RowMajMat::Constant(row, col, 1E-8)); + velocities.setEpsilon(1E-8); + + Advection advection = Advection(concentrations, velocities); + + advection.setPorosity(RowMajMat::Constant(row, col, 0.2)); + advection.setIterations(1); + advection.setTimestep(10000); // create boundaries - Boundary bcH = Boundary(gridHeads); - bcH.setBoundarySideConstant(BC_SIDE_LEFT, 10); + Boundary &bcH = velocities.getBoundaryConditions(); + bcH.setBoundarySideConstant(BC_SIDE_LEFT, 1); bcH.setBoundarySideConstant(BC_SIDE_RIGHT, 1); - bcH.setBoundarySideClosed(BC_SIDE_TOP); - bcH.setBoundarySideClosed(BC_SIDE_BOTTOM); - Boundary bcC = Boundary(gridConc); - bcC.setBoundarySideConstant(BC_SIDE_LEFT, 0.1); - bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 1); + bcH.setBoundarySideConstant(BC_SIDE_TOP, 1); + bcH.setBoundarySideConstant(BC_SIDE_BOTTOM, 1); + // bcH.setBoundarySideClosed(BC_SIDE_TOP); + // bcH.setBoundarySideClosed(BC_SIDE_BOTTOM); + bcH.setInnerBoundary(row / 2, col / 2, 10); + // + Boundary &bcC = advection.getBoundaryConditions(); + // bcC.setBoundarySideConstant(BC_SIDE_LEFT, 0.1); + // bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 1); bcC.setBoundarySideClosed(BC_SIDE_TOP); bcC.setBoundarySideClosed(BC_SIDE_BOTTOM); - Velocities velocities = Velocities(gridHeads, bcH); - velocities.setOutputCSV(CSV_OUTPUT_ON); - velocities.setK(1E-2); - velocities.setEpsilon(1E-8); - velocities.setInjh(10); - velocities.setIterations(0); - // calculate steady hydraulic heads - velocities.run(); - - std::cout << "Velocities simulation finished." << std::endl; - std::cout << hydHeads << std::endl; - - // set true for steady case - Advection advection = Advection(velocities, gridConc, bcC, true); - advection.setPorosity(0.2); - advection.setIterations(21); - // set timestep close almost exactly to clf to test advection - advection.setTimestep(5039.05); - // velocities.setOutputCSV(CSV_OUTPUT_VERBOSE); advection.run(); + std::cout << velocities.getConcentrationMatrix() << std::endl; + std::cout << velocities.getVelocitiesX() << std::endl; + std::cout << velocities.getVelocitiesY() << std::endl; + std::cout << "Advection simulation finished." << std::endl; std::cout << concentrations << std::endl; -} \ No newline at end of file +} diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 8924b9c..c6f2eb8 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,9 +1,9 @@ -add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp) -add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp) -add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp) +# add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp) +# add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp) +# add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp) add_executable(Advection Advection.cpp) -target_link_libraries(BTCS_2D_proto_example tug) -target_link_libraries(FTCS_2D_proto_closed_mdl tug) -target_link_libraries(FTCS_2D_proto_example_mdl tug) +# target_link_libraries(BTCS_2D_proto_example tug) +# target_link_libraries(FTCS_2D_proto_closed_mdl tug) +# target_link_libraries(FTCS_2D_proto_example_mdl tug) target_link_libraries(Advection tug) diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp index 855e1b5..ff448c5 100644 --- a/include/tug/Advection/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -8,18 +8,16 @@ #pragma once #include "tug/Core/Matrix.hpp" +#include +#include #include -#include -#include -#include +#include #include #include #include #include -#include #include -#include #include #include #include @@ -29,7 +27,15 @@ using namespace Eigen; namespace tug { -template class Advection { +template +class Advection : public BaseSimulationGrid { +private: + T timestep{-1}; + int numThreads{omp_get_num_procs()}; + + Velocities &velocities; + RowMajMat porosity; + public: /** * @brief Construct a new Advection object, used to calculate material @@ -46,53 +52,45 @@ public: * @param Steady Used to choose between Steady and Transient case. Either true * or false */ - Advection(Velocities &_velocities, Grid &_grid, Boundary &_bc, - bool Steady) - : velocities(_velocities), grid(_grid), bc(_bc), - outx(_velocities.getOutx()), outy(_velocities.getOuty()), - Steady(Steady) {}; + Advection(RowMajMat &origin, + Velocities &_velocities) + : BaseSimulationGrid(origin), velocities(_velocities) { + tug_assert(origin.rows() == velocities.getConcentrationMatrix().rows() && + origin.cols() == velocities.getConcentrationMatrix().cols(), + "Advection grid and Velocities must have the same dimensions"); + }; + + Advection(T *data, std::size_t rows, std::size_t cols, + Velocities &_velocities) + : BaseSimulationGrid(data, rows, cols), velocities(_velocities) { + tug_assert(rows == velocities.getConcentrationMatrix().rows() && + cols == velocities.getConcentrationMatrix().cols(), + "Advection grid and Velocities must have the same dimensions"); + }; /** * @brief Sets the porosity of the medium * * @param porosity new porosity value */ - void setPorosity(T porosity) { - if (porosity < 0 || porosity > 1) { - throw std::invalid_argument( - "Porosity must be a value between 0 and 1 (inclusive)"); - } + void setPorosity(const RowMajMat &porosity) { + tug_assert(porosity.rows() == this->rows() && + porosity.cols() == this->cols(), + "Porosity matrix must have the same dimensions as the grid"); + tug_assert(porosity.minCoeff() >= 0 && porosity.maxCoeff() <= 1, + "Porosity must be a value between 0 and 1 (inclusive)"); + this->porosity = porosity; } - /** - * @brief Set the desired iterations to be calculated. A value greater - * than zero must be specified here. Setting iterations is required. - * - * @param iterations Number of iterations to be simulated. - */ - void setIterations(int iterations) { - if (iterations <= 0) { - throw std::invalid_argument( - "Number of iterations must be greater than zero. Provided value: " + - std::to_string(iterations)); - } - this->iterations = iterations; - }; - /** * @brief Set the size of the timestep. Must be greater than zero * * @param timestep Size of the timestep */ void setTimestep(T timestep) { - if (timestep <= 0) { - throw std::invalid_argument( - "Timestep must be greater than zero. Provided value: " + - std::to_string(timestep)); - } else { - this->timestep = timestep; - } + tug_assert(timestep > 0, "Timestep must be greater than zero"); + this->timestep = timestep; } /** @@ -118,113 +116,54 @@ public: } } - /** - * @brief Set the option to output the results to a CSV file. Off by default. - * - * - * @param csv_output Valid output option. The following options can be set - * here: - * - CSV_OUTPUT_OFF: do not produce csv output - * - CSV_OUTPUT_ON: produce csv output with last - * concentration matrix - if (csv_output < CSV_OUTPUT_OFF || csv_output > CSV_OUTPUT_VERBOSE) { - throw std::invalid_argument("Invalid CSV output option given: " + - std::to_string(csv_output)); - } - void setOutputCSV(CSV_OUTPUT csv_output) { - if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) { - throw std::invalid_argument("Invalid CSV output option given!"); - } + void run() { + this->setDomain(velocities.domainX(), velocities.domainY()); - this->csv_output = csv_output; - } - - /** - * @brief Creates a CSV file with a name containing the current simulation - * parameters. If the data name already exists, an additional counter - * is appended to the name. The name of the file is built up as follows: - * + + + - * +.csv - * - * @return string Filename with configured simulation parameters. - */ - std::string createCSVfile(std::string Type) const { - std::ofstream file; - int appendIdent = 0; - std::string appendIdentString; - - std::string row = std::to_string(grid.getRow()); - std::string col = std::to_string(grid.getCol()); - std::string numIterations = std::to_string(iterations); - - std::string filename = - Type + "_" + row + "_" + col + "_" + numIterations + ".csv"; - - while (std::filesystem::exists(filename)) { - appendIdent += 1; - appendIdentString = std::to_string(appendIdent); - filename = Type + "_" + row + "_" + col + "_" + numIterations + "-" + - appendIdentString + ".csv"; - } - - file.open(filename); - if (!file) { - throw std::runtime_error("Failed to open file: " + filename); - } - - file.close(); - - return filename; - } - /** - * @brief Writes the currently calculated Concentration values of the grid - * into the CSV file with the passed filename. - * - * @param filename Name of the file to which the Concentration values are - * to be written. - */ - void printConcentrationCSV(const std::string &filename) const { - std::ofstream file; - - file.open(filename, std::ios_base::app); - if (!file) { - throw std::runtime_error("Failed to open file: " + filename); - } - - Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols); - file << grid.getConcentrations().format(do_not_align) << std::endl; - file << std::endl << std::endl; - file.close(); + velocities.run(); + + for (int i = 0; i < this->getIterations(); i++) { + if constexpr (hyd_mode == HYDRAULIC_MODE::TRANSIENT) { + velocities.run(); + } + adv(); + } } +private: /** * @brief function calculating material transport for one timestep */ void adv() { - int rows = grid.getRow(); - int cols = grid.getCol(); - T volume = grid.getDeltaRow() * grid.getDeltaCol(); + const std::size_t rows = this->rows(); + const std::size_t cols = this->cols(); + const T volume = this->deltaCol() * this->deltaRow(); - RowMajMat &newConcentrations = grid.getConcentrations(); + RowMajMatMap &newConcentrations = this->getConcentrationMatrix(); + + const RowMajMat &outx = this->velocities.getVelocitiesX(); + const RowMajMat &outy = this->velocities.getVelocitiesY(); + const Boundary &bc = this->getBoundaryConditions(); // Calculate Courant-Levy-Frederich condition - T maxFx = std::max(abs(outx.maxCoeff()), abs(outx.minCoeff())); - T maxFy = std::max(abs(outy.maxCoeff()), abs(outy.minCoeff())); - T maxF = std::max(maxFx, maxFy); + const T maxFx = std::max(abs(outx.maxCoeff()), abs(outx.minCoeff())); + const T maxFy = std::max(abs(outy.maxCoeff()), abs(outy.minCoeff())); + const T maxF = std::max(maxFx, maxFy); - if (maxF == 0) { - throw std::runtime_error("Division by zero: maxF is zero."); - } + tug_assert(maxF != 0, "Division by zero: maxF is zero."); - T cvf = abs((volume * porosity) / maxF); - int innerSteps = (int)ceil(timestep / cvf); - T innerTimestep = timestep / innerSteps; + const RowMajMat volumeTimesPorosity = volume * porosity; + + const T cvf = (volumeTimesPorosity / maxF).maxCoeff(); + const int innerSteps = (int)ceil(timestep / cvf); + const T innerTimestep = timestep / innerSteps; + + const RowMajMat multiplier = volumeTimesPorosity * (1 / innerTimestep); for (int k = 0; k < innerSteps; k++) { - const Eigen::MatrixX oldConcentrations = newConcentrations; + const RowMajMat oldConcentrations = newConcentrations; // Calculate sum of incoming/outgoing Flow*Concentration in x-direction in each // cell -#pragma omp parallel for num_threads(numThreads) schedule(static) +#pragma omp parallel for num_threads(numThreads) for (int i = 0; i < rows; i++) { for (int j = 0; j < cols + 1; j++) { if (j == 0) { @@ -318,51 +257,15 @@ public: } } } - newConcentrations = - oldConcentrations + - newConcentrations * (innerTimestep / (porosity * volume)); + + for (std::size_t row_i = 0; row_i < rows; row_i++) { + for (std::size_t col_i = 0; col_i < cols; col_i++) { + newConcentrations(row_i, col_i) = + oldConcentrations(row_i, col_i) + + newConcentrations(row_i, col_i) * multiplier(row_i, col_i); + } + } } } - - void run() { - std::string filename; - if (csv_output >= CSV_OUTPUT_ON) { - filename = createCSVfile("Concentrations"); - } - - if (Steady == false) { - velocities.setTimestep(timestep); - velocities.setIterations(1); - } - - for (int i = 0; i < iterations; i++) { - if (csv_output >= CSV_OUTPUT_VERBOSE) { - printConcentrationCSV(filename); - } - // if steady==false update charge and velocities with equal timestep - if (Steady == false) { - velocities.run(); - } - adv(); - } - - if (csv_output >= CSV_OUTPUT_ON) { - printConcentrationCSV(filename); - } - } - -private: - Grid &grid; - Boundary &bc; - Velocities &velocities; - bool Steady{true}; - int iterations{-1}; - int innerIterations{1}; - T timestep{-1}; - int numThreads{omp_get_num_procs()}; - T porosity{1}; - const Eigen::MatrixX &outx; - const Eigen::MatrixX &outy; - CSV_OUTPUT csv_output{CSV_OUTPUT_OFF}; }; } // namespace tug diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp index bcb78ba..0c50671 100644 --- a/include/tug/Core/Numeric/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -82,14 +82,13 @@ template static void FTCS_1D(SimulationInput &input) { const T ×tep = input.timestep; RowMajMatMap &concentrations_grid = input.concentrations; - // matrix for concentrations at time t+1 - RowMajMat concentrations_t1 = concentrations_grid; const auto &alphaX = input.alphaX; const auto &bc = input.boundaries; - checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax, input.colMax); + // matrix for concentrations at time t+1 + RowMajMat concentrations_t1 = concentrations_grid; // only one row in 1D case -> row constant at index 0 int row = 0; @@ -99,6 +98,9 @@ template static void FTCS_1D(SimulationInput &input) { // inner cells // independent of boundary condition type for (int col = 1; col < colMax - 1; col++) { + if (inner_bc[col].first) { + continue; + } const T &conc_c = concentrations_grid(row, col); const T &conc_left = concentrations_grid(row, col - 1); const T &conc_right = concentrations_grid(row, col + 1); @@ -115,8 +117,8 @@ template static void FTCS_1D(SimulationInput &input) { } // left boundary; hold column constant at index 0 - { - int col = 0; + int col = 0; + if (!inner_bc[col].first) { const T &conc_c = concentrations_grid(row, col); const T &conc_right = concentrations_grid(row, col + 1); const T &alpha_c = alphaX(row, col); @@ -132,8 +134,8 @@ template static void FTCS_1D(SimulationInput &input) { } // right boundary; hold column constant at max index - { - int col = colMax - 1; + col = colMax - 1; + if (!inner_bc[col].first) { const T &conc_c = concentrations_grid(row, col); const T &conc_left = concentrations_grid(row, col - 1); const T &alpha_c = alphaX(row, col); @@ -161,27 +163,28 @@ static void FTCS_2D(SimulationInput &input, int numThreads) { const T ×tep = input.timestep; RowMajMatMap &concentrations_grid = input.concentrations; - - // matrix for concentrations at time t+1 - RowMajMat concentrations_t1 = concentrations_grid; - const auto &alphaX = input.alphaX; const auto &alphaY = input.alphaY; const auto &bc = input.boundaries; + const T sx = timestep / (deltaCol * deltaCol); + const T sy = timestep / (deltaRow * deltaRow); + checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax, input.colMax); - const T sx = timestep / (deltaCol * deltaCol); - const T sy = timestep / (deltaRow * deltaRow); + // matrix for concentrations at time t+1 + RowMajMat concentrations_t1 = concentrations_grid; #pragma omp parallel for num_threads(numThreads) for (std::size_t row_i = 0; row_i < rowMax; row_i++) { for (std::size_t col_i = 0; col_i < colMax; col_i++) { + if (bc.getInnerBoundary(row_i, col_i).first) { + continue; + } // horizontal change T horizontal_change; { - const T &conc_c = concentrations_grid(row_i, col_i); const T &alpha_c = alphaX(row_i, col_i); diff --git a/include/tug/tug.hpp b/include/tug/tug.hpp index bc797ae..6770de3 100644 --- a/include/tug/tug.hpp +++ b/include/tug/tug.hpp @@ -1,6 +1,7 @@ #pragma once -// #include +#include +#include #include #include #include From d8c8a734aa373ba7e3b6117fc8288a76f4a6acc8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 6 Feb 2025 14:55:54 +0100 Subject: [PATCH 10/20] test(diffusion): Verify symmetry in diffusion simulation --- test/testDiffusion.cpp | 44 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 75a51b7..b1faf7c 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -214,3 +214,47 @@ DIFFUSION_TEST(ConstantInnerCell) { EXPECT_FALSE((concentrations_result.array() < 0.0).any()); } + +DIFFUSION_TEST(Symmetry) { + // Arrange + constexpr std::size_t rows = 25; + constexpr std::size_t cols = 25; + + constexpr std::size_t center_row = rows / 2; + constexpr std::size_t center_col = cols / 2; + + tug::RowMajMat concentrations = + tug::RowMajMat::Constant(rows, cols, 1); + + tug::RowMajMat alpha = + tug::RowMajMat::Constant(rows, cols, 1E-2); + + tug::Diffusion sim(concentrations); + + sim.setDomain(100, 100); + sim.setAlphaX(alpha); + sim.setAlphaY(alpha); + // choose a high number of iterations, which lead to small changes in ULP + // between symmetric cells + sim.setIterations(10000); + sim.setTimestep(10); + + tug::Boundary &bcH = sim.getBoundaryConditions(); + + bcH.setInnerBoundary(center_row, center_col, 10); + + sim.run(); + + // check symmetry + for (std::size_t i_rows = 0; i_rows <= center_row; i_rows++) { + for (std::size_t i_cols = 0; i_cols <= center_col; i_cols++) { + if (i_rows == center_row && i_cols == center_col) { + continue; + } + // to avoid floating point errors, we check with ASSERT_DOUBLE_EQ with a + // precision of ULP(4), see https://stackoverflow.com/a/4149599 + ASSERT_DOUBLE_EQ(concentrations(i_rows, i_cols), + concentrations(rows - i_rows - 1, cols - i_cols - 1)); + } + } +} From 16b361c85b4f263d1c7459839a17bf5789e9634c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 6 Feb 2025 16:18:04 +0100 Subject: [PATCH 11/20] fix(velocities): prevent redundant velocity calculations --- include/tug/Advection/Advection.hpp | 10 +++++++--- include/tug/Advection/Velocities.hpp | 6 ++++++ test/testVelocities.cpp | 12 +++--------- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp index ff448c5..6fd52b7 100644 --- a/include/tug/Advection/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -119,13 +119,17 @@ public: void run() { this->setDomain(velocities.domainX(), velocities.domainY()); - velocities.run(); + if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) { + if (!velocities.isSteady()) { + velocities.run(); + } + } for (int i = 0; i < this->getIterations(); i++) { if constexpr (hyd_mode == HYDRAULIC_MODE::TRANSIENT) { velocities.run(); } - adv(); + startAdvection(); } } @@ -133,7 +137,7 @@ private: /** * @brief function calculating material transport for one timestep */ - void adv() { + void startAdvection() { const std::size_t rows = this->rows(); const std::size_t cols = this->cols(); const T volume = this->deltaCol() * this->deltaRow(); diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp index 8cf2fc3..b98a195 100644 --- a/include/tug/Advection/Velocities.hpp +++ b/include/tug/Advection/Velocities.hpp @@ -43,6 +43,7 @@ private: T timestep{-1}; T epsilon{1E-5}; int numThreads{omp_get_num_procs()}; + bool steady{false}; RowMajMat velocitiesX; RowMajMat velocitiesY; @@ -118,6 +119,8 @@ public: this->permKY = alphaY; } + bool isSteady() const { return steady; } + /** * @brief Set the timestep per iteration * @@ -219,6 +222,9 @@ public: oldConcentrations = this->getConcentrationMatrix(); (void)calculate_hydraulic_flow(input); } while (!checkConvergance(oldConcentrations)); + + steady = true; + } else { if (timestep == -1) { throw_invalid_argument("Timestep is not set"); diff --git a/test/testVelocities.cpp b/test/testVelocities.cpp index 0c8af66..40f16d2 100644 --- a/test/testVelocities.cpp +++ b/test/testVelocities.cpp @@ -20,9 +20,7 @@ VELOCITIES_TEST(SteadyStateCenter) { tug::RowMajMat hydHeads = tug::RowMajMat::Constant(rows, cols, 1); - tug::RowMajMat permKX = - tug::RowMajMat::Constant(rows, cols, K); - tug::RowMajMat permKY = + tug::RowMajMat permK = tug::RowMajMat::Constant(rows, cols, K); tug::Velocities &bcH = velo.getBoundaryConditions(); - bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1); - bcH.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1); - bcH.setBoundarySideConstant(tug::BC_SIDE_TOP, 1); - bcH.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1); bcH.setInnerBoundary(center_row, center_col, 10); From 2250ee3f6f297936dcae03b8c258d0d2d678dc95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 7 Feb 2025 09:51:46 +0100 Subject: [PATCH 12/20] refactor(advection): move steady state check to velocities --- include/tug/Advection/Advection.hpp | 4 +--- include/tug/Advection/Velocities.hpp | 4 ++++ 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp index 6fd52b7..a20de13 100644 --- a/include/tug/Advection/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -120,9 +120,7 @@ public: this->setDomain(velocities.domainX(), velocities.domainY()); if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) { - if (!velocities.isSteady()) { - velocities.run(); - } + velocities.run(); } for (int i = 0; i < this->getIterations(); i++) { diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp index b98a195..f62f8d9 100644 --- a/include/tug/Advection/Velocities.hpp +++ b/include/tug/Advection/Velocities.hpp @@ -206,6 +206,10 @@ public: .deltaCol = this->deltaCol()}; if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) { + if (steady) { + return; + } + const T deltaColSquare = this->deltaCol() * this->deltaCol(); const T deltaRowSquare = this->deltaRow() * this->deltaRow(); const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); From bdb44b4fd55fc5294e18e784d5be8aa1a8d9ba4b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 7 Feb 2025 09:52:43 +0100 Subject: [PATCH 13/20] fix(ftcs): add return statement for undefined boundary condition --- include/tug/Core/Numeric/FTCS.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp index 0c50671..28d2db6 100644 --- a/include/tug/Core/Numeric/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -54,6 +54,7 @@ constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center, } tug_assert(false, "Undefined Boundary Condition Type!"); + return 0; } template @@ -71,8 +72,6 @@ checkAndSetConstantInnerCells(const Boundary &bc, const auto &col = rowcol.second; concentrations(row, col) = value; } - - return; } // FTCS solution for 1D grid From 031905b4c88e6c20a725fd9e448868ecaa4f055e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 7 Feb 2025 09:53:00 +0100 Subject: [PATCH 14/20] test: Add advection test case with left-to-right flow --- test/CMakeLists.txt | 3 +- test/testAdvection.cpp | 69 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 71 insertions(+), 1 deletion(-) create mode 100644 test/testAdvection.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a932356..c0d2a51 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -12,9 +12,10 @@ FetchContent_MakeAvailable(googletest) add_executable(testTug setup.cpp testDiffusion.cpp +testVelocities.cpp +testAdvection.cpp testFTCS.cpp testBoundary.cpp -testVelocities.cpp ) target_link_libraries(testTug tug GTest::gtest) diff --git a/test/testAdvection.cpp b/test/testAdvection.cpp new file mode 100644 index 0000000..ac4c2db --- /dev/null +++ b/test/testAdvection.cpp @@ -0,0 +1,69 @@ +#include +#include + +#define ADVECTION_TEST(x) TEST(Advection, x) + +ADVECTION_TEST(LeftToRight) { + constexpr std::size_t rows = 21; + constexpr std::size_t cols = 21; + + constexpr double K = 1E-2; + constexpr double timestep = 5039.05; + constexpr std::size_t iterations = 21; + constexpr double porosity = 0.2; + + constexpr double epsilon = 1E-13; + + tug::RowMajMat hydHeads = + tug::RowMajMat::Constant(rows, cols, 1); + tug::RowMajMat concentrations = + tug::RowMajMat::Constant(rows, cols, 0); + + tug::RowMajMat permK = + tug::RowMajMat::Constant(rows, cols, K); + + tug::Velocities + velocities(hydHeads); + velocities.setDomain(100, 100); + velocities.setPermKX(permK); + velocities.setPermKY(permK); + velocities.setEpsilon(1E-8); + + tug::Advection advection(concentrations, velocities); + + advection.setPorosity(tug::RowMajMat::Constant(rows, cols, porosity)); + advection.setIterations(iterations); + advection.setTimestep(timestep); + + tug::Boundary &bcH = velocities.getBoundaryConditions(); + bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 10); + bcH.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1); + + tug::Boundary &bcC = advection.getBoundaryConditions(); + bcC.setBoundarySideConstant(tug::BC_SIDE_LEFT, 0.1); + bcC.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1); + + advection.run(); + + // check if the concentration is transported from left to right + for (std::size_t i_rows = 0; i_rows < rows; i_rows++) { + for (std::size_t i_cols = 0; i_cols < cols - 1; i_cols++) { + if (i_cols == 0) { + EXPECT_LE(concentrations(i_rows, i_cols), 10); + } else { + EXPECT_GE(concentrations(i_rows, i_cols), + concentrations(i_rows, i_cols + 1)); + } + } + } + + // the values should also be equal from top to bottom + for (std::size_t i_cols = 0; i_cols < cols; i_cols++) { + const double &ref = concentrations(0, i_cols); + for (std::size_t i_rows = 1; i_rows < rows; i_rows++) { + // check if the values are equal within the epsilon range + EXPECT_NEAR(ref, concentrations(i_rows, i_cols), epsilon); + } + } +} From 2be7b82f70894f2b54411ab3b44f8e23e9481e40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 7 Feb 2025 13:24:13 +0100 Subject: [PATCH 15/20] feat: Apply inner boundary conditions before simulation steps --- include/tug/Advection/Advection.hpp | 2 ++ include/tug/Advection/Velocities.hpp | 2 ++ include/tug/Core/BaseSimulation.hpp | 12 ++++++++++++ include/tug/Core/Numeric/FTCS.hpp | 22 ---------------------- include/tug/Diffusion/Diffusion.hpp | 2 ++ 5 files changed, 18 insertions(+), 22 deletions(-) diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp index a20de13..8409279 100644 --- a/include/tug/Advection/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -119,6 +119,8 @@ public: void run() { this->setDomain(velocities.domainX(), velocities.domainY()); + this->applyInnerBoundaries(); + if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) { velocities.run(); } diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp index f62f8d9..fd9580f 100644 --- a/include/tug/Advection/Velocities.hpp +++ b/include/tug/Advection/Velocities.hpp @@ -193,6 +193,8 @@ public: // if iterations < 1 calculate hydraulic charge until steady state is // reached + this->applyInnerBoundaries(); + SimulationInput input = {.concentrations = this->getConcentrationMatrix(), .alphaX = this->getPermKX(), diff --git a/include/tug/Core/BaseSimulation.hpp b/include/tug/Core/BaseSimulation.hpp index 82a1480..3dc8268 100644 --- a/include/tug/Core/BaseSimulation.hpp +++ b/include/tug/Core/BaseSimulation.hpp @@ -64,6 +64,18 @@ private: T delta_col; T delta_row; +protected: + void applyInnerBoundaries() { + const auto &inner_bc = boundaryConditions.getInnerBoundaries(); + if (inner_bc.empty()) { + return; + } + + for (const auto &[rowcol, value] : inner_bc) { + concentrationMatrix(rowcol.first, rowcol.second) = value; + } + } + public: /** * @brief Constructs a BaseSimulationGrid from a given RowMajMat object. diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp index 28d2db6..4ccf1a4 100644 --- a/include/tug/Core/Numeric/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -57,23 +57,6 @@ constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center, return 0; } -template -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; - } - - for (const auto &[rowcol, value] : inner_boundaries) { - const auto &row = rowcol.first; - const auto &col = rowcol.second; - concentrations(row, col) = value; - } -} - // FTCS solution for 1D grid template static void FTCS_1D(SimulationInput &input) { const std::size_t &colMax = input.colMax; @@ -84,8 +67,6 @@ 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); // matrix for concentrations at time t+1 RowMajMat concentrations_t1 = concentrations_grid; @@ -169,9 +150,6 @@ static void FTCS_2D(SimulationInput &input, int numThreads) { const T sx = timestep / (deltaCol * deltaCol); const T sy = timestep / (deltaRow * deltaRow); - checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax, - input.colMax); - // matrix for concentrations at time t+1 RowMajMat concentrations_t1 = concentrations_grid; diff --git a/include/tug/Diffusion/Diffusion.hpp b/include/tug/Diffusion/Diffusion.hpp index e353550..f977514 100644 --- a/include/tug/Diffusion/Diffusion.hpp +++ b/include/tug/Diffusion/Diffusion.hpp @@ -347,6 +347,8 @@ public: auto begin = std::chrono::high_resolution_clock::now(); + this->applyInnerBoundaries(); + SimulationInput sim_input = {.concentrations = this->getConcentrationMatrix(), .alphaX = this->getAlphaX(), From 8b273a59b16bcef7b07cc6f1b02f3d16854f0e80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 7 Feb 2025 14:38:26 +0100 Subject: [PATCH 16/20] [wip] --- examples/Advection.cpp | 60 +++-- include/tug/Advection/Advection.hpp | 338 +++++++++++++++++++-------- include/tug/Advection/Velocities.hpp | 3 + 3 files changed, 267 insertions(+), 134 deletions(-) diff --git a/examples/Advection.cpp b/examples/Advection.cpp index c2ca66a..8a6d192 100644 --- a/examples/Advection.cpp +++ b/examples/Advection.cpp @@ -1,7 +1,5 @@ -#include "tug/Advection/Advection.hpp" -#include "tug/Advection/Velocities.hpp" -#include "tug/Core/Matrix.hpp" -#include +#include "tug/Boundary.hpp" +#include #include #include @@ -9,53 +7,53 @@ using namespace Eigen; using namespace tug; int main(int argc, char *argv[]) { - int row = 5; - int col = 5; - // create two grids of equal size, grid1 for hydraulics heads, grid2 for - // Concentrations + constexpr std::size_t row = 5; + constexpr std::size_t col = 5; RowMajMat hydHeads = RowMajMat::Constant(row, col, 1); RowMajMat concentrations = RowMajMat::Constant(row, col, 0); - hydHeads(row / 2, col / 2) = 10; - concentrations(row / 2, col / 2) = 1; Velocities velocities(hydHeads); - velocities.setDomain(100, 100); - velocities.setPermKX(RowMajMat::Constant(row, col, 1E-8)); - velocities.setPermKY(RowMajMat::Constant(row, col, 1E-8)); + + velocities.setDomain(1, 1); + velocities.setPermKX(RowMajMat::Constant(row, col, 3E-7)); + velocities.setPermKY(RowMajMat::Constant(row, col, 3E-7)); velocities.setEpsilon(1E-8); Advection advection = Advection(concentrations, velocities); advection.setPorosity(RowMajMat::Constant(row, col, 0.2)); - advection.setIterations(1); - advection.setTimestep(10000); + advection.setIterations(6); + // 1 hour + advection.setTimestep(6666); // create boundaries Boundary &bcH = velocities.getBoundaryConditions(); - bcH.setBoundarySideConstant(BC_SIDE_LEFT, 1); - bcH.setBoundarySideConstant(BC_SIDE_RIGHT, 1); - bcH.setBoundarySideConstant(BC_SIDE_TOP, 1); - bcH.setBoundarySideConstant(BC_SIDE_BOTTOM, 1); - // bcH.setBoundarySideClosed(BC_SIDE_TOP); - // bcH.setBoundarySideClosed(BC_SIDE_BOTTOM); - bcH.setInnerBoundary(row / 2, col / 2, 10); - // + bcH.setBoundarySideConstant(BC_SIDE_LEFT, 10); + bcH.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + // bcH.setBoundarySideConstant(BC_SIDE_TOP, 1); + // bcH.setBoundarySideConstant(BC_SIDE_BOTTOM, 1); + // bcH.setInnerBoundary(row / 2, col / 2, 10); + Boundary &bcC = advection.getBoundaryConditions(); - // bcC.setBoundarySideConstant(BC_SIDE_LEFT, 0.1); - // bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 1); - bcC.setBoundarySideClosed(BC_SIDE_TOP); - bcC.setBoundarySideClosed(BC_SIDE_BOTTOM); + bcC.setBoundarySideConstant(BC_SIDE_LEFT, 1); + bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + // bcC.setInnerBoundary(row / 2, col / 2, 1); + // bcC.setBoundarySideConstant(BC_SIDE_LEFT, 0); + // bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + // bcC.setBoundarySideConstant(BC_SIDE_TOP, 0); + // bcC.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); advection.run(); - std::cout << velocities.getConcentrationMatrix() << std::endl; - std::cout << velocities.getVelocitiesX() << std::endl; - std::cout << velocities.getVelocitiesY() << std::endl; + std::cout << velocities.getConcentrationMatrix() << std::endl << std::endl; - std::cout << "Advection simulation finished." << std::endl; + std::cout << velocities.getVelocitiesX() << std::endl + << std::endl + << velocities.getVelocitiesY() << std::endl + << std::endl; std::cout << concentrations << std::endl; } diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp index 8409279..bb9a572 100644 --- a/include/tug/Advection/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -12,6 +12,8 @@ #include #include #include +#include +#include #include #include #include @@ -140,135 +142,265 @@ private: void startAdvection() { 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 T volume = this->deltaCol() * this->deltaRow(); RowMajMatMap &newConcentrations = this->getConcentrationMatrix(); - const RowMajMat &outx = this->velocities.getVelocitiesX(); - const RowMajMat &outy = this->velocities.getVelocitiesY(); + RowMajMat veloX = this->velocities.getVelocitiesX(); + RowMajMat veloY = this->velocities.getVelocitiesY(); const Boundary &bc = this->getBoundaryConditions(); // Calculate Courant-Levy-Frederich condition - const T maxFx = std::max(abs(outx.maxCoeff()), abs(outx.minCoeff())); - const T maxFy = std::max(abs(outy.maxCoeff()), abs(outy.minCoeff())); + const T maxFx = std::max(abs(veloX.maxCoeff()), abs(veloX.minCoeff())); + const T maxFy = std::max(abs(veloY.maxCoeff()), abs(veloY.minCoeff())); const T maxF = std::max(maxFx, maxFy); tug_assert(maxF != 0, "Division by zero: maxF is zero."); - const RowMajMat volumeTimesPorosity = volume * porosity; + // TODO: sustitute 1 by porisity + const RowMajMat porevolume = + RowMajMat::Constant(rows, cols, volume) * 1; - const T cvf = (volumeTimesPorosity / maxF).maxCoeff(); - const int innerSteps = (int)ceil(timestep / cvf); + const T cfl = (porevolume / maxF).maxCoeff() / 2; + const int innerSteps = (int)ceil(timestep / cfl); const T innerTimestep = timestep / innerSteps; - const RowMajMat multiplier = volumeTimesPorosity * (1 / innerTimestep); + std::cout << "CFL (adv) timestep: " << cfl << std::endl; + std::cout << "Inner iterations (adv): " << innerSteps << std::endl; - for (int k = 0; k < innerSteps; k++) { + // const RowMajMat multiplier = porevolume * (1 / innerTimestep); + + // convienient renaming to same memory + auto &fluxX = veloX; + auto &fluxY = veloY; + + for (std::size_t k = 0; k < innerSteps; k++) { const RowMajMat oldConcentrations = newConcentrations; -// Calculate sum of incoming/outgoing Flow*Concentration in x-direction in each -// cell -#pragma omp parallel for num_threads(numThreads) - for (int i = 0; i < rows; i++) { - for (int j = 0; j < cols + 1; j++) { - if (j == 0) { - if (bc.getBoundaryElementType(BC_SIDE_LEFT, i) != BC_TYPE_CLOSED) { - if (outx(i, j) > 0) { - // outx positive -> flow from border to cell i,j - newConcentrations(i, j) += - outx(i, j) * bc.getBoundaryElementValue(BC_SIDE_LEFT, i); - } else if (outx(i, j) < 0) { - // outx negative -> flow from i,j towards border - newConcentrations(i, j) += outx(i, j) * oldConcentrations(i, j); - } - } - } else if (j == cols) { - if (bc.getBoundaryElementType(BC_SIDE_RIGHT, i) != BC_TYPE_CLOSED) { - if (outx(i, j) > 0) { - // outx positive-> flow from i,j-1 towards border - newConcentrations(i, j - 1) -= - outx(i, j) * oldConcentrations(i, j - 1); - } else if (outx(i, j) < 0) { - // outx negative -> flow from border to cell i,j-1 - newConcentrations(i, j - 1) -= - outx(i, j) * bc.getBoundaryElementValue(BC_SIDE_LEFT, i); - } - } - } - // flow between inner cells - else { - // outx positive -> flow from cell i,j-1 towards cell i,j - if (outx(i, j) > 0) { - newConcentrations(i, j - 1) -= - outx(i, j) * oldConcentrations(i, j - 1); - newConcentrations(i, j) += - outx(i, j) * oldConcentrations(i, j - 1); - } - // outx negative -> flow from cell i,j toward cell i,j-1 - else if (outx(i, j) < 0) { - newConcentrations(i, j - 1) -= - outx(i, j) * oldConcentrations(i, j); - newConcentrations(i, j) += outx(i, j) * oldConcentrations(i, j); - } - } - } - } -// calculate sum in y-direction -// parallelize outer loop over columns to ensure thread-safety, each thread only -// modifies cells within its column -#pragma omp parallel for num_threads(numThreads) - for (int j = 0; j < cols; j++) { - for (int i = 0; i < rows + 1; i++) { - if (i == 0) { - if (bc.getBoundaryElementType(BC_SIDE_TOP, j) != BC_TYPE_CLOSED) { - if (outy(i, j) > 0) { - // outy positive -> flow from border to cell i,j - newConcentrations(i, j) += - outy(i, j) * bc.getBoundaryElementValue(BC_SIDE_TOP, j); - } else if (outy(i, j) < 0) { - // outy negative -> flow from i,j towards border - newConcentrations(i, j) += outy(i, j) * oldConcentrations(i, j); - } - } - } else if (i == rows) { - if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, j) != - BC_TYPE_CLOSED) { - if (outy(i, j) > 0) { - // outy positive-> flow from i-1,j towards border - newConcentrations(i - 1, j) -= - outy(i, j) * oldConcentrations(i - 1, j); - } else if (outy(i, j) < 0) { - // outy negative -> flow from border to cell i,j-1 - newConcentrations(i - 1, j) -= - outy(i, j) * bc.getBoundaryElementValue(BC_SIDE_BOTTOM, j); - } - } - } - // flow between inner cells - else { - // outy positive -> flow from cell i-1,j towards cell i,j - if (outy(i, j) > 0) { - newConcentrations(i - 1, j) -= - outy(i, j) * oldConcentrations(i - 1, j); - newConcentrations(i, j) += - outy(i, j) * oldConcentrations(i - 1, j); - } - // outy negative -> flow from cell i,j toward cell i-1,j - else if (outy(i, j) < 0) { - newConcentrations(i - 1, j) -= - outy(i, j) * oldConcentrations(i, j); - newConcentrations(i, j) += outy(i, j) * oldConcentrations(i, j); + + // Update flux with concentrations in x-direction + for (std::size_t row_i = 0; row_i < fluxX.rows(); row_i++) { + for (std::size_t col_i = 0; col_i < fluxX.cols(); col_i++) { + T ¤tFlux = fluxX(row_i, col_i); + const std::int32_t cellIndex = + (fluxX(row_i, col_i) > 0) ? col_i - 1 : col_i; + + if (cellIndex < 0 || cellIndex >= cols) { + const auto bcElement = bc.getBoundaryElement( + cellIndex < 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i); + // if we have a constant boundary, we use the value of the boundary + if (bcElement.getType() == BC_TYPE_CONSTANT) { + currentFlux *= bcElement.getValue() * (2 * innerTimestep) / dx; + continue; } + + // otherwise it is a closed boundary + currentFlux = 0; + continue; } + // otherwise we use the concentration of the inflow/outflow cell, + // multiplied by the velocity + currentFlux *= + oldConcentrations(row_i, cellIndex) * innerTimestep / dx; } } + // Update flux with concentrations in y-direction + // for (std::size_t row_i = 0; row_i < fluxY.rows(); row_i++) { + // for (std::size_t col_i = 0; col_i < fluxY.cols(); col_i++) { + // T ¤tFlux = fluxY(row_i, col_i); + // const std::int32_t cellIndex = + // (fluxY(row_i, col_i) > 0) ? row_i - 1 : row_i; + // + // if (cellIndex < 0 || cellIndex >= rows) { + // const auto bcElement = bc.getBoundaryElement( + // cellIndex < 0 ? BC_SIDE_TOP : BC_SIDE_BOTTOM, col_i); + // // if we have a constant boundary, we use the value of the + // boundary if (bcElement.getType() == BC_TYPE_CONSTANT) { + // currentFlux *= bcElement.getValue() * (2 * innerTimestep) / dy; + // continue; + // } + // + // // otherwise it is a closed boundary + // currentFlux = 0; + // continue; + // } + // // otherwise we use the concentration of the inflow/outflow cell, + // // multiplied by the velocity + // currentFlux *= + // oldConcentrations(cellIndex, col_i) * innerTimestep / dy; + // } + // } + + // Update concentrations for (std::size_t row_i = 0; row_i < rows; row_i++) { for (std::size_t col_i = 0; col_i < cols; col_i++) { + const T horizontalFlux = + fluxX(row_i, col_i) - fluxX(row_i, col_i + 1); + const T verticalFlux = 0; newConcentrations(row_i, col_i) = - oldConcentrations(row_i, col_i) + - newConcentrations(row_i, col_i) * multiplier(row_i, col_i); + oldConcentrations(row_i, col_i) + horizontalFlux + verticalFlux; } } + + // const Boundary leftCorner = bc.getBoundaryElement(BC_SIDE_LEFT, 0); + // newConcentrations(0, 0) = oldConcentrations(0, 0) + ; + // + // for (std::size_t col_i = 0; col_i < cols; col_i++) { + // T horizontalFlux = oldConcentrations(0, col_i); + // if (col_i == 0) { + // switch (bc.getBoundaryElementType(BC_SIDE_LEFT, 0)) { + // case BC_TYPE_CLOSED: { + // break; + // } + // case BC_TYPE_CONSTANT: { + // newConcentrations(0, 0) = + // oldConcentrations(0, 0) + + // (bc.getBoundaryElementValue(BC_SIDE_LEFT, 0) * veloX(0, 0)) * + // (2 * innerTimestep) / dx; + // newConcentrations(rows - 1, 0) = + // oldConcentrations(rows - 1, 0) + + // (bc.getBoundaryElementValue(BC_SIDE_LEFT, 0) * + // veloX(rows - 1, 0)) * + // (2 * innerTimestep) / dx; + // break; + // } + // } + // newConcentrations(0, 0) -= + // veloX(0, 1) * oldConcentrations(0, 1) * innerTimestep / dx; + // } + // } + // + // for (std::size_t row_i = 1; row_i < rows - 1; row_i++) { + // for (std::size_t col_i = 1; col_i < cols - 1; col_i++) { + // const T fluxHorizontal = + // veloX(row_i, col_i) * oldConcentrations(row_i, col_i - 1) - + // veloX(row_i, col_i + 1) * oldConcentrations(row_i, col_i); + // const T fluxVertical = + // veloY(row_i, col_i) * oldConcentrations(row_i - 1, col_i) - + // veloY(row_i + 1, col_i) * oldConcentrations(row_i, col_i); + // newConcentrations(row_i, col_i) = + // oldConcentrations(row_i, col_i) + + // (fluxHorizontal / dx + fluxVertical / dy) * innerTimestep; + // } + // } + + // Calculate sum of incoming/outgoing Flow*Concentration in x-direction in + // each cell #pragma omp parallel for num_threads(numThreads) + // for (int i = 0; i < rows; i++) { + // for (int j = 0; j < cols + 1; j++) { + // if (j == 0) { + // if (bc.getBoundaryElementType(BC_SIDE_LEFT, i) != + // BC_TYPE_CLOSED) { + // if (veloX(i, j) > 0) { + // // outx positive -> flow from border to cell i,j + // newConcentrations(i, j) += + // veloX(i, j) * + // bc.getBoundaryElementValue(BC_SIDE_LEFT, i); + // } else if (veloX(i, j) < 0) { + // // outx negative -> flow from i,j towards border + // newConcentrations(i, j) += veloX(i, j) * + // oldConcentrations(i, j); + // } + // } + // } else if (j == cols) { + // if (bc.getBoundaryElementType(BC_SIDE_RIGHT, i) != + // BC_TYPE_CLOSED) { + // if (veloX(i, j) > 0) { + // // outx positive-> flow from i,j-1 towards border + // newConcentrations(i, j - 1) -= + // veloX(i, j) * oldConcentrations(i, j - 1); + // } else if (veloX(i, j) < 0) { + // // outx negative -> flow from border to cell i,j-1 + // newConcentrations(i, j - 1) -= + // veloX(i, j) * + // bc.getBoundaryElementValue(BC_SIDE_LEFT, i); + // } + // } + // } + // // flow between inner cells + // else { + // // outx positive -> flow from cell i,j-1 towards cell i,j + // if (veloX(i, j) > 0) { + // newConcentrations(i, j - 1) -= + // veloX(i, j) * oldConcentrations(i, j - 1); + // newConcentrations(i, j) += + // veloX(i, j) * oldConcentrations(i, j - 1); + // } + // // outx negative -> flow from cell i,j toward cell i,j-1 + // else if (veloX(i, j) < 0) { + // newConcentrations(i, j - 1) += + // veloX(i, j) * oldConcentrations(i, j); + // newConcentrations(i, j) += veloX(i, j) * + // oldConcentrations(i, j); + // } + // } + // } + // } + // // calculate sum in y-direction + // // parallelize outer loop over columns to ensure thread-safety, each + // thread only + // // modifies cells within its column + // #pragma omp parallel for num_threads(numThreads) + // for (int j = 0; j < cols; j++) { + // for (int i = 0; i < rows + 1; i++) { + // if (i == 0) { + // if (bc.getBoundaryElementType(BC_SIDE_TOP, j) != + // BC_TYPE_CLOSED) { + // if (veloY(i, j) > 0) { + // // outy positive -> flow from border to cell i,j + // newConcentrations(i, j) += + // veloY(i, j) * + // bc.getBoundaryElementValue(BC_SIDE_TOP, j); + // } else if (veloY(i, j) < 0) { + // // outy negative -> flow from i,j towards border + // newConcentrations(i, j) += veloY(i, j) * + // oldConcentrations(i, j); + // } + // } + // } else if (i == rows) { + // if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, j) != + // BC_TYPE_CLOSED) { + // if (veloY(i, j) > 0) { + // // outy positive-> flow from i-1,j towards border + // newConcentrations(i - 1, j) -= + // veloY(i, j) * oldConcentrations(i - 1, j); + // } else if (veloY(i, j) < 0) { + // // outy negative -> flow from border to cell i,j-1 + // newConcentrations(i - 1, j) -= + // veloY(i, j) * + // bc.getBoundaryElementValue(BC_SIDE_BOTTOM, j); + // } + // } + // } + // // flow between inner cells + // else { + // // outy positive -> flow from cell i-1,j towards cell i,j + // if (veloY(i, j) > 0) { + // newConcentrations(i - 1, j) -= + // veloY(i, j) * oldConcentrations(i - 1, j); + // newConcentrations(i, j) += + // veloY(i, j) * oldConcentrations(i - 1, j); + // } + // // outy negative -> flow from cell i,j toward cell i-1,j + // else if (veloY(i, j) < 0) { + // newConcentrations(i - 1, j) -= + // veloY(i, j) * oldConcentrations(i, j); + // newConcentrations(i, j) += veloY(i, j) * + // oldConcentrations(i, j); + // } + // } + // } + // } + + // for (std::size_t row_i = 0; row_i < rows; row_i++) { + // for (std::size_t col_i = 0; col_i < cols; col_i++) { + // newConcentrations(row_i, col_i) = + // oldConcentrations(row_i, col_i) + + // newConcentrations(row_i, col_i) * multiplier(row_i, col_i); + // } + // } } } }; diff --git a/include/tug/Advection/Velocities.hpp b/include/tug/Advection/Velocities.hpp index fd9580f..6679db8 100644 --- a/include/tug/Advection/Velocities.hpp +++ b/include/tug/Advection/Velocities.hpp @@ -219,6 +219,9 @@ public: const T maxK = std::max(this->permKX.maxCoeff(), this->permKY.maxCoeff()); // Calculate largest possible timestep, depending on K and gridsize + std::cout << "CFL (hydHead) timestep: " << minDeltaSquare / (4 * maxK) + << std::endl; + setTimestep(minDeltaSquare / (4 * maxK)); input.timestep = this->timestep; From 1391716ecb2556e579b225c190435b4001b62586 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Fri, 7 Feb 2025 17:28:52 +0100 Subject: [PATCH 17/20] [wip] fix advection scheme --- examples/Advection.cpp | 6 +++--- include/tug/Advection/Advection.hpp | 9 +++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/examples/Advection.cpp b/examples/Advection.cpp index 8a6d192..2ce374f 100644 --- a/examples/Advection.cpp +++ b/examples/Advection.cpp @@ -17,7 +17,7 @@ int main(int argc, char *argv[]) { tug::HYDRAULIC_RESOLVE::EXPLICIT> velocities(hydHeads); - velocities.setDomain(1, 1); + velocities.setDomain(5, 5); velocities.setPermKX(RowMajMat::Constant(row, col, 3E-7)); velocities.setPermKY(RowMajMat::Constant(row, col, 3E-7)); velocities.setEpsilon(1E-8); @@ -25,9 +25,9 @@ int main(int argc, char *argv[]) { Advection advection = Advection(concentrations, velocities); advection.setPorosity(RowMajMat::Constant(row, col, 0.2)); - advection.setIterations(6); + advection.setIterations(3); // 1 hour - advection.setTimestep(6666); + advection.setTimestep(1.6666e+06); // create boundaries Boundary &bcH = velocities.getBoundaryConditions(); diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp index bb9a572..a73d9be 100644 --- a/include/tug/Advection/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -155,20 +155,21 @@ private: // Calculate Courant-Levy-Frederich condition const T maxFx = std::max(abs(veloX.maxCoeff()), abs(veloX.minCoeff())); const T maxFy = std::max(abs(veloY.maxCoeff()), abs(veloY.minCoeff())); - const T maxF = std::max(maxFx, maxFy); + const T maxFlux = std::max(maxFx, maxFy); - tug_assert(maxF != 0, "Division by zero: maxF is zero."); + tug_assert(maxFlux != 0, "Division by zero: maxF is zero."); // TODO: sustitute 1 by porisity const RowMajMat porevolume = RowMajMat::Constant(rows, cols, volume) * 1; - const T cfl = (porevolume / maxF).maxCoeff() / 2; + const T cfl = (porevolume / maxFlux).minCoeff(); const int innerSteps = (int)ceil(timestep / cfl); const T innerTimestep = timestep / innerSteps; std::cout << "CFL (adv) timestep: " << cfl << std::endl; std::cout << "Inner iterations (adv): " << innerSteps << std::endl; + std::cout << "Inner timestep (adv): " << innerTimestep << std::endl; // const RowMajMat multiplier = porevolume * (1 / innerTimestep); @@ -191,7 +192,7 @@ private: cellIndex < 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i); // if we have a constant boundary, we use the value of the boundary if (bcElement.getType() == BC_TYPE_CONSTANT) { - currentFlux *= bcElement.getValue() * (2 * innerTimestep) / dx; + currentFlux *= bcElement.getValue() * (innerTimestep) / dx; continue; } From 1ce20c972c7f7d2b1d889bce864661c7a019d5c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 11 Feb 2025 07:50:32 +0100 Subject: [PATCH 18/20] fix(advection): correct flux calculation with velocities If two or more inner iterations were required, instead of velocities the previous calculated flux was used as velocity. This lead to erroneous results. --- include/tug/Advection/Advection.hpp | 69 +++++++++++++++-------------- 1 file changed, 36 insertions(+), 33 deletions(-) diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp index a73d9be..7970275 100644 --- a/include/tug/Advection/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -13,7 +13,6 @@ #include #include #include -#include #include #include #include @@ -173,9 +172,8 @@ private: // const RowMajMat multiplier = porevolume * (1 / innerTimestep); - // convienient renaming to same memory - auto &fluxX = veloX; - auto &fluxY = veloY; + RowMajMat fluxX(rows, cols + 1); + RowMajMat fluxY(rows + 1, cols); for (std::size_t k = 0; k < innerSteps; k++) { const RowMajMat oldConcentrations = newConcentrations; @@ -185,14 +183,16 @@ private: for (std::size_t col_i = 0; col_i < fluxX.cols(); col_i++) { T ¤tFlux = fluxX(row_i, col_i); const std::int32_t cellIndex = - (fluxX(row_i, col_i) > 0) ? col_i - 1 : col_i; + (veloX(row_i, col_i) > 0) ? col_i - 1 : col_i; + const T ¤tVelo = veloX(row_i, col_i); if (cellIndex < 0 || cellIndex >= cols) { const auto bcElement = bc.getBoundaryElement( cellIndex < 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i); // if we have a constant boundary, we use the value of the boundary if (bcElement.getType() == BC_TYPE_CONSTANT) { - currentFlux *= bcElement.getValue() * (innerTimestep) / dx; + currentFlux = + currentVelo * bcElement.getValue() * (innerTimestep) / dx; continue; } @@ -202,37 +202,40 @@ private: } // otherwise we use the concentration of the inflow/outflow cell, // multiplied by the velocity - currentFlux *= - oldConcentrations(row_i, cellIndex) * innerTimestep / dx; + currentFlux = currentVelo * oldConcentrations(row_i, cellIndex) * + innerTimestep / dx; } } // Update flux with concentrations in y-direction - // for (std::size_t row_i = 0; row_i < fluxY.rows(); row_i++) { - // for (std::size_t col_i = 0; col_i < fluxY.cols(); col_i++) { - // T ¤tFlux = fluxY(row_i, col_i); - // const std::int32_t cellIndex = - // (fluxY(row_i, col_i) > 0) ? row_i - 1 : row_i; - // - // if (cellIndex < 0 || cellIndex >= rows) { - // const auto bcElement = bc.getBoundaryElement( - // cellIndex < 0 ? BC_SIDE_TOP : BC_SIDE_BOTTOM, col_i); - // // if we have a constant boundary, we use the value of the - // boundary if (bcElement.getType() == BC_TYPE_CONSTANT) { - // currentFlux *= bcElement.getValue() * (2 * innerTimestep) / dy; - // continue; - // } - // - // // otherwise it is a closed boundary - // currentFlux = 0; - // continue; - // } - // // otherwise we use the concentration of the inflow/outflow cell, - // // multiplied by the velocity - // currentFlux *= - // oldConcentrations(cellIndex, col_i) * innerTimestep / dy; - // } - // } + for (std::size_t row_i = 0; row_i < fluxY.rows(); row_i++) { + for (std::size_t col_i = 0; col_i < fluxY.cols(); col_i++) { + T ¤tFlux = fluxY(row_i, col_i); + const std::int32_t cellIndex = + (veloY(row_i, col_i) > 0) ? row_i - 1 : row_i; + const T ¤tVelo = veloY(row_i, col_i); + + if (cellIndex < 0 || cellIndex >= rows) { + const auto bcElement = bc.getBoundaryElement( + cellIndex < 0 ? BC_SIDE_TOP : BC_SIDE_BOTTOM, col_i); + // if we have a constant boundary, we use the value of the + // boundary + if (bcElement.getType() == BC_TYPE_CONSTANT) { + currentFlux = + currentVelo * bcElement.getValue() * (2 * innerTimestep) / dy; + continue; + } + + // otherwise it is a closed boundary + currentFlux = 0; + continue; + } + // otherwise we use the concentration of the inflow/outflow cell, + // multiplied by the velocity + currentFlux = currentVelo * oldConcentrations(cellIndex, col_i) * + innerTimestep / dy; + } + } // Update concentrations for (std::size_t row_i = 0; row_i < rows; row_i++) { From da8973674e3b4dc8cd6733748889d96af6494564 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 11 Feb 2025 11:01:48 +0100 Subject: [PATCH 19/20] fix: Correct flux calculation and boundary condition handling --- include/tug/Advection/Advection.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp index 7970275..88563a9 100644 --- a/include/tug/Advection/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -182,9 +182,8 @@ private: for (std::size_t row_i = 0; row_i < fluxX.rows(); row_i++) { for (std::size_t col_i = 0; col_i < fluxX.cols(); col_i++) { T ¤tFlux = fluxX(row_i, col_i); - const std::int32_t cellIndex = - (veloX(row_i, col_i) > 0) ? col_i - 1 : col_i; const T ¤tVelo = veloX(row_i, col_i); + const std::int32_t cellIndex = (currentVelo > 0) ? col_i - 1 : col_i; if (cellIndex < 0 || cellIndex >= cols) { const auto bcElement = bc.getBoundaryElement( @@ -211,9 +210,8 @@ private: for (std::size_t row_i = 0; row_i < fluxY.rows(); row_i++) { for (std::size_t col_i = 0; col_i < fluxY.cols(); col_i++) { T ¤tFlux = fluxY(row_i, col_i); - const std::int32_t cellIndex = - (veloY(row_i, col_i) > 0) ? row_i - 1 : row_i; const T ¤tVelo = veloY(row_i, col_i); + const std::int32_t cellIndex = (currentVelo > 0) ? row_i - 1 : row_i; if (cellIndex < 0 || cellIndex >= rows) { const auto bcElement = bc.getBoundaryElement( @@ -242,7 +240,7 @@ private: for (std::size_t col_i = 0; col_i < cols; col_i++) { const T horizontalFlux = fluxX(row_i, col_i) - fluxX(row_i, col_i + 1); - const T verticalFlux = 0; + const T verticalFlux = fluxY(row_i, col_i) - fluxY(row_i + 1, col_i); newConcentrations(row_i, col_i) = oldConcentrations(row_i, col_i) + horizontalFlux + verticalFlux; } From f2f4d6fca8bd8bcd0ad84c3a15cdff6cbab387dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 11 Feb 2025 14:17:06 +0100 Subject: [PATCH 20/20] fix(advection): correct flux calculation in advection scheme --- include/tug/Advection/Advection.hpp | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/include/tug/Advection/Advection.hpp b/include/tug/Advection/Advection.hpp index 88563a9..58ed9b1 100644 --- a/include/tug/Advection/Advection.hpp +++ b/include/tug/Advection/Advection.hpp @@ -9,6 +9,7 @@ #include "tug/Core/Matrix.hpp" #include +#include #include #include #include @@ -190,8 +191,7 @@ private: cellIndex < 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i); // if we have a constant boundary, we use the value of the boundary if (bcElement.getType() == BC_TYPE_CONSTANT) { - currentFlux = - currentVelo * bcElement.getValue() * (innerTimestep) / dx; + currentFlux = currentVelo * bcElement.getValue(); continue; } @@ -201,8 +201,7 @@ private: } // otherwise we use the concentration of the inflow/outflow cell, // multiplied by the velocity - currentFlux = currentVelo * oldConcentrations(row_i, cellIndex) * - innerTimestep / dx; + currentFlux = currentVelo * oldConcentrations(row_i, cellIndex); } } @@ -219,8 +218,7 @@ private: // if we have a constant boundary, we use the value of the // boundary if (bcElement.getType() == BC_TYPE_CONSTANT) { - currentFlux = - currentVelo * bcElement.getValue() * (2 * innerTimestep) / dy; + currentFlux = currentVelo * bcElement.getValue(); continue; } @@ -230,8 +228,7 @@ private: } // otherwise we use the concentration of the inflow/outflow cell, // multiplied by the velocity - currentFlux = currentVelo * oldConcentrations(cellIndex, col_i) * - innerTimestep / dy; + currentFlux = currentVelo * oldConcentrations(cellIndex, col_i); } } @@ -241,8 +238,10 @@ private: const T horizontalFlux = fluxX(row_i, col_i) - fluxX(row_i, col_i + 1); const T verticalFlux = fluxY(row_i, col_i) - fluxY(row_i + 1, col_i); - newConcentrations(row_i, col_i) = - oldConcentrations(row_i, col_i) + horizontalFlux + verticalFlux; + newConcentrations(row_i, col_i) = oldConcentrations(row_i, col_i) + + innerTimestep / + porevolume(row_i, col_i) * + (horizontalFlux + verticalFlux); } }