max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-16 02:48:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

feat: Implement steady-state hydraulic charge calculation

Browse Source
This commit is contained in:
Max Lübke 2025-02-05 15:42:36 +01:00
parent ca94cebba2
commit 7a1d9bb5b7
6 changed files with 133 additions and 108 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

160
include/tug/Advection/Velocities.hpp
View File

@ -8,6 +8,7 @@
#pragma once #pragma once
#include "tug/Core/Numeric/SimulationInput.hpp"
#include <algorithm> #include <algorithm>
#include <cstddef> #include <cstddef>
#include <iostream> #include <iostream>
@ -16,12 +17,12 @@
#include <string> #include <string>
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
#include <tug/Core/BaseSimulation.hpp>
#include <tug/Core/Matrix.hpp> #include <tug/Core/Matrix.hpp>
#include <tug/Core/Numeric/BTCS.hpp> #include <tug/Core/Numeric/BTCS.hpp>
#include <tug/Core/Numeric/FTCS.hpp> #include <tug/Core/Numeric/FTCS.hpp>
#include <tug/Core/TugUtils.hpp> #include <tug/Core/TugUtils.hpp>
#include <tug/Diffusion/Diffusion.hpp> #include <tug/Diffusion/Diffusion.hpp>
#include <tug/Grid.hpp>
#ifdef _OPENMP #ifdef _OPENMP
#include <omp.h> #include <omp.h>
@ -32,38 +33,38 @@
using namespace Eigen; using namespace Eigen;
namespace tug { namespace tug {
enum HYDRAULIC_MODE { TRANSIENT, STEADY_STATE }; enum class HYDRAULIC_MODE { TRANSIENT, STEADY_STATE };
enum HYDRAULIC_RESOLVE { EXPLICIT, IMPLICIT }; enum class HYDRAULIC_RESOLVE { EXPLICIT, IMPLICIT };
template <class T, HYDRAULIC_MODE hyd_mode, HYDRAULIC_RESOLVE hyd_resolve> template <class T, HYDRAULIC_MODE hyd_mode, HYDRAULIC_RESOLVE hyd_resolve>
class Velocities : public BaseSimulation { class Velocities : public BaseSimulationGrid<T> {
private: private:
int innerIterations{1}; int innerIterations{1};
T timestep{-1}; T timestep{-1};
T epsilon{1E-5}; T epsilon{1E-5};
int numThreads{omp_get_num_procs()}; int numThreads{omp_get_num_procs()};
Grid<T> &grid;
Boundary<T> &bc;
RowMajMat<T> velocitiesX; RowMajMat<T> velocitiesX;
RowMajMat<T> velocitiesY; RowMajMat<T> velocitiesY;
RowMajMat<T> alphaX;
RowMajMat<T> alphaY;
public: public:
/** Velocities(RowMajMat<T> &origin)
* @brief Construct a new Velocities object, used to calculate Hydraulic : BaseSimulationGrid<T>(origin),
* Charge and Darcy-Velocities. A timestep and a number of iterations can be velocitiesX(origin.rows(), origin.cols() + 1),
* set. By default iterations is set to -1. If the number of iterations is set velocitiesY(origin.rows() + 1, origin.cols()),
* to a value below 1 the simulation will run until the Hydraulic Charge alphaX(origin.rows(), origin.cols()),
* converges. The Epsilon value to check convergence can be set, the default alphaY(origin.rows(), origin.cols()) {};
* is 1E-5. CSV Output is off by default.
* Velocities(T *data, std::size_t rows, std::size_t cols)
* @param grid Valid grid object : BaseSimulationGrid<T>(data, rows, cols), velocitiesX(rows, cols + 1),
* @param bc Valid boundary condition object velocitiesY(rows + 1, cols), alphaX(rows, cols), alphaY(rows, cols) {};
*/
Velocities(Grid<T> &_grid, Boundary<T> &_bc) Velocities(T *data, std::size_t length)
: grid(_grid), bc(_bc), velocitiesX(grid.getRow(), grid.getCol() + 1), : BaseSimulationGrid<T>(data, 1, length), velocitiesX(1, length + 1),
velocitiesY(grid.getRow() + 1, grid.getCol()) {}; alphaX(1, length) {};
/** /**
* @brief Set the epsilon value, the relativ error allowed for convergence * @brief Set the epsilon value, the relativ error allowed for convergence
@ -77,22 +78,61 @@ public:
this->epsilon = epsilon; this->epsilon = epsilon;
} }
/**
* @brief Get the alphaX matrix.
*
* @return RowMajMat<T>& Reference to the alphaX matrix.
*/
RowMajMat<T> &getAlphaX() { return alphaX; }
/**
* @brief Get the alphaY matrix.
*
* @return RowMajMat<T>& Reference to the alphaY matrix.
*/
RowMajMat<T> &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<T> &alphaX) { this->alphaX = alphaX; }
/**
* @brief Set the alphaY matrix.
*
* @param alphaY The new alphaY matrix.
*/
void setAlphaY(const RowMajMat<T> &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 * @brief Set the timestep per iteration
* *
* @param timestep timestep per iteration * @param timestep timestep per iteration
*/ */
void setTimestep(T timestep) { void setTimestep(T timestep) override {
if (timestep <= 0) { if (timestep <= 0) {
throw std::invalid_argument("Timestep must be greater than zero"); throw std::invalid_argument("Timestep must be greater than zero");
} }
this->timestep = timestep; this->timestep = timestep;
const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); const T deltaColSquare = this->deltaCol() * this->deltaCol();
const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); const T deltaRowSquare = this->deltaRow() * this->deltaRow();
const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare); const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
const T maxK = const T maxK = std::max(this->alphaX.maxCoeff(), this->alphaY.maxCoeff());
std::max(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff());
T cfl = minDeltaSquare / (4 * maxK); T cfl = minDeltaSquare / (4 * maxK);
if (timestep > cfl) { if (timestep > cfl) {
@ -146,32 +186,45 @@ public:
* @brief Simulation of hydraulic charge either until convergence, * @brief Simulation of hydraulic charge either until convergence,
* or for a number of selected timesteps. Calculation of Darcy-velocities. * 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 // if iterations < 1 calculate hydraulic charge until steady state is
// reached // reached
if constexpr (hyd_mode == STEADY_STATE) { SimulationInput<T> input = {.concentrations =
// Calculate largest possible timestep, depending on K and gridsize this->getConcentrationMatrix(),
const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol(); .alphaX = this->getAlphaX(),
const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow(); .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 minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
const T maxK = const T maxK = std::max(this->alphaX.maxCoeff(), this->alphaY.maxCoeff());
std::max(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff()); // Calculate largest possible timestep, depending on K and gridsize
setTimestep(minDeltaSquare / (4 * maxK)); setTimestep(minDeltaSquare / (4 * maxK));
input.timestep = this->timestep;
RowMajMat<T> oldConcentrations; RowMajMat<T> oldConcentrations;
do { do {
oldConcentrations = grid.getConcentrations(); oldConcentrations = this->getConcentrationMatrix();
(void)calculate_hydraulic_flow(); (void)calculate_hydraulic_flow(input);
} while (!checkConvergance(oldConcentrations, grid.getConcentrations())); } while (!checkConvergance(oldConcentrations));
} else { } else {
if (timestep == -1) { if (timestep == -1) {
throw_invalid_argument("Timestep is not set"); throw_invalid_argument("Timestep is not set");
} }
for (int i = 0; i < innerIterations; i++) { 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 * @brief Calculate the new hydraulic charge using FTCS
*/ */
void calculate_hydraulic_flow() { void calculate_hydraulic_flow(SimulationInput<T> &sim_in) {
if constexpr (hyd_resolve == EXPLICIT) { if constexpr (hyd_resolve == HYDRAULIC_RESOLVE::EXPLICIT) {
FTCS_2D(this->grid, this->bc, this->timestep, this->numThreads); FTCS_2D(sim_in, numThreads);
} else { } 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 * containing old and new Charge values, the relative error is below the
* selected Epsilon * selected Epsilon
*/ */
bool checkConvergance(const RowMajMat<T> &oldHeads, bool checkConvergance(const RowMajMat<T> &oldHeads) {
const RowMajMat<T> &newHeads) { const auto abs_err = (oldHeads - this->getConcentrationMatrix()).cwiseAbs();
const auto abs_err = (oldHeads - newHeads).cwiseAbs(); const auto rel_err = abs_err.cwiseQuotient(this->getConcentrationMatrix());
const auto rel_err = abs_err.cwiseQuotient(newHeads);
return rel_err.maxCoeff() < epsilon; return rel_err.maxCoeff() < epsilon;
} }
@ -210,14 +262,16 @@ private:
* directions * directions
*/ */
void computeFluidVelocities() { void computeFluidVelocities() {
const std::size_t rows = grid.getRow(); const std::size_t rows = this->rows();
const std::size_t cols = grid.getCol(); const std::size_t cols = this->cols();
const T dx = grid.getDeltaRow(); const T dx = this->deltaCol();
const T dy = grid.getDeltaCol(); const T dy = this->deltaRow();
const RowMajMat<T> &hydraulicCharges = grid.getConcentrations(); const RowMajMat<T> &hydraulicCharges = this->getConcentrationMatrix();
const RowMajMat<T> &permKX = grid.getAlphaX(); const RowMajMat<T> &permKX = this->alphaX;
const RowMajMat<T> &permKY = grid.getAlphaY(); const RowMajMat<T> &permKY = this->alphaY;
const Boundary<T> &bc = this->getBoundaryConditions();
// calculate velocities in x-direction // calculate velocities in x-direction
for (std::size_t i_rows = 0; i_rows < rows; i_rows++) { for (std::size_t i_rows = 0; i_rows < rows; i_rows++) {
@ -229,7 +283,7 @@ private:
} }
case BC_TYPE_CONSTANT: { case BC_TYPE_CONSTANT: {
velocitiesX(i_rows, 0) = velocitiesX(i_rows, 0) =
-permKX(i_rows, 0) * -this->alphaX(i_rows, 0) *
(hydraulicCharges(i_rows, 0) - bc_left.getValue()) / (dx / 2); (hydraulicCharges(i_rows, 0) - bc_left.getValue()) / (dx / 2);
break; break;
} }

20
include/tug/Core/Numeric/FTCS.hpp
View File

@ -8,13 +8,13 @@
#ifndef FTCS_H_ #ifndef FTCS_H_
#define FTCS_H_ #define FTCS_H_
#include "tug/Core/TugUtils.hpp" #include "tug/Core/Matrix.hpp"
#include <cstddef> #include <cstddef>
#include <cstring> #include <cstring>
#include <iterator>
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
#include <tug/Core/Matrix.hpp> #include <tug/Core/Matrix.hpp>
#include <tug/Core/Numeric/SimulationInput.hpp> #include <tug/Core/Numeric/SimulationInput.hpp>
#include <tug/Core/TugUtils.hpp>
#ifdef _OPENMP #ifdef _OPENMP
#include <omp.h> #include <omp.h>
@ -57,17 +57,15 @@ constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center,
} }
template <typename T> template <typename T>
static inline void checkAndSetConstantInnerCells(const Boundary<T> &bc, static inline void
Grid<T> &grid) { checkAndSetConstantInnerCells(const Boundary<T> &bc,
RowMajMatMap<T> &concentrations, std::size_t rows,
std::size_t cols) {
const auto &inner_boundaries = bc.getInnerBoundaries(); const auto &inner_boundaries = bc.getInnerBoundaries();
if (inner_boundaries.empty()) { if (inner_boundaries.empty()) {
return; return;
} }
auto &concentrations = grid.getConcentrations();
const auto rows = grid.getRow();
const auto cols = grid.getCol();
for (const auto &[rowcol, value] : inner_boundaries) { for (const auto &[rowcol, value] : inner_boundaries) {
const auto &row = rowcol.first; const auto &row = rowcol.first;
const auto &col = rowcol.second; const auto &col = rowcol.second;
@ -90,6 +88,9 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
const auto &alphaX = input.alphaX; const auto &alphaX = input.alphaX;
const auto &bc = input.boundaries; const auto &bc = input.boundaries;
checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax,
input.colMax);
// only one row in 1D case -> row constant at index 0 // only one row in 1D case -> row constant at index 0
int row = 0; int row = 0;
@ -168,6 +169,9 @@ static void FTCS_2D(SimulationInput<T> &input, int numThreads) {
const auto &alphaY = input.alphaY; const auto &alphaY = input.alphaY;
const auto &bc = input.boundaries; const auto &bc = input.boundaries;
checkAndSetConstantInnerCells(bc, concentrations_grid, input.rowMax,
input.colMax);
const T sx = timestep / (deltaCol * deltaCol); const T sx = timestep / (deltaCol * deltaCol);
const T sy = timestep / (deltaRow * deltaRow); const T sy = timestep / (deltaRow * deltaRow);

3
include/tug/Core/Numeric/SimulationInput.hpp
View File

@ -11,9 +11,8 @@ template <typename T> struct SimulationInput {
const RowMajMat<T> &alphaX; const RowMajMat<T> &alphaX;
const RowMajMat<T> &alphaY; const RowMajMat<T> &alphaY;
const Boundary<T> boundaries; const Boundary<T> boundaries;
const std::uint8_t dim; const std::uint8_t dim;
const T timestep; T timestep;
const std::size_t rowMax; const std::size_t rowMax;
const std::size_t colMax; const std::size_t colMax;
const T deltaRow; const T deltaRow;

1
include/tug/tug.hpp
View File

@ -4,4 +4,3 @@
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
#include <tug/Core/Matrix.hpp> #include <tug/Core/Matrix.hpp>
#include <tug/Diffusion/Diffusion.hpp> #include <tug/Diffusion/Diffusion.hpp>
#include <tug/Grid.hpp>

33
test/testDiffusion.cpp
View File

@ -3,8 +3,7 @@
#include "gtest/gtest.h" #include "gtest/gtest.h"
#include <gtest/gtest.h> #include <gtest/gtest.h>
#include <stdexcept> #include <stdexcept>
#include <tug/Diffusion/Diffusion.hpp>
#include <tug/tug.hpp>
#include <Eigen/src/Core/Matrix.h> #include <Eigen/src/Core/Matrix.h>
#include <string> #include <string>
@ -215,33 +214,3 @@ DIFFUSION_TEST(ConstantInnerCell) {
EXPECT_FALSE((concentrations_result.array() < 0.0).any()); 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<double, tug::FTCS_APPROACH> 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());
}

24
test/testVelocities.cpp
View File

@ -25,11 +25,15 @@ VELOCITIES_TEST(SteadyStateCenter) {
tug::RowMajMat<double> permKY = tug::RowMajMat<double> permKY =
tug::RowMajMat<double>::Constant(rows, cols, K); tug::RowMajMat<double>::Constant(rows, cols, K);
tug::Grid64 gridHeads(hydHeads); tug::Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
gridHeads.setDomain(100, 100); tug::HYDRAULIC_RESOLVE::EXPLICIT>
gridHeads.setAlpha(permKX, permKY); velo(hydHeads);
tug::Boundary bcH = tug::Boundary(gridHeads); velo.setDomain(100, 100);
velo.setAlphaX(permKX);
velo.setAlphaY(permKY);
tug::Boundary<double> &bcH = velo.getBoundaryConditions();
bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1); bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1);
bcH.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1); bcH.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1);
bcH.setBoundarySideConstant(tug::BC_SIDE_TOP, 1); bcH.setBoundarySideConstant(tug::BC_SIDE_TOP, 1);
@ -37,14 +41,10 @@ VELOCITIES_TEST(SteadyStateCenter) {
bcH.setInnerBoundary(center_row, center_col, 10); bcH.setInnerBoundary(center_row, center_col, 10);
tug::Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE, velo.run();
tug::HYDRAULIC_RESOLVE::EXPLICIT>
velocities(gridHeads, bcH);
velocities.run(); const auto &velocitiesX = velo.getVelocitiesX();
const auto &velocitiesY = velo.getVelocitiesY();
const auto &velocitiesX = velocities.getVelocitiesX();
const auto &velocitiesY = velocities.getVelocitiesY();
// Assert // Assert
@ -81,4 +81,4 @@ VELOCITIES_TEST(SteadyStateCenter) {
} }
} }
} }
} }