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

[wip]

Browse Source
This commit is contained in:
Max Lübke 2025-02-07 14:38:26 +01:00
parent 2be7b82f70
commit 8b273a59b1
3 changed files with 267 additions and 134 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
examples/Advection.cpp
View File

@ -1,7 +1,5 @@
#include "tug/Advection/Advection.hpp"
#include "tug/Advection/Velocities.hpp"
#include "tug/Core/Matrix.hpp"
#include <Eigen/Eigen>
#include "tug/Boundary.hpp"
#include <cstddef>
#include <iostream>
#include <tug/tug.hpp>
@ -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<double> hydHeads = RowMajMat<double>::Constant(row, col, 1);
RowMajMat<double> concentrations = RowMajMat<double>::Constant(row, col, 0);
hydHeads(row / 2, col / 2) = 10;
concentrations(row / 2, col / 2) = 1;
Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
tug::HYDRAULIC_RESOLVE::EXPLICIT>
velocities(hydHeads);
velocities.setDomain(100, 100);
velocities.setPermKX(RowMajMat<double>::Constant(row, col, 1E-8));
velocities.setPermKY(RowMajMat<double>::Constant(row, col, 1E-8));
velocities.setDomain(1, 1);
velocities.setPermKX(RowMajMat<double>::Constant(row, col, 3E-7));
velocities.setPermKY(RowMajMat<double>::Constant(row, col, 3E-7));
velocities.setEpsilon(1E-8);
Advection advection = Advection(concentrations, velocities);
advection.setPorosity(RowMajMat<double>::Constant(row, col, 0.2));
advection.setIterations(1);
advection.setTimestep(10000);
advection.setIterations(6);
// 1 hour
advection.setTimestep(6666);
// create boundaries
Boundary<double> &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<double> &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;
}

338
include/tug/Advection/Advection.hpp
View File

@ -12,6 +12,8 @@
#include <Eigen/src/Core/Matrix.h>
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <iterator>
#include <stdexcept>
#include <stdlib.h>
#include <string>
@ -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<T> &newConcentrations = this->getConcentrationMatrix();
const RowMajMat<T> &outx = this->velocities.getVelocitiesX();
const RowMajMat<T> &outy = this->velocities.getVelocitiesY();
RowMajMat<T> veloX = this->velocities.getVelocitiesX();
RowMajMat<T> veloY = this->velocities.getVelocitiesY();
const Boundary<T> &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<T> volumeTimesPorosity = volume * porosity;
// TODO: sustitute 1 by porisity
const RowMajMat<T> porevolume =
RowMajMat<T>::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<T> 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<T> 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<T> 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 &currentFlux = 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 &currentFlux = 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<T> 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);
// }
// }
}
}
};

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

@ -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;