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] [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;