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