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

fix(advection): correct flux calculation in advection scheme

Browse Source
This commit is contained in:
Max Lübke 2025-02-11 14:17:06 +01:00
parent da8973674e
commit f2f4d6fca8
1 changed files with 9 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -9,6 +9,7 @@
#include "tug/Core/Matrix.hpp" #include "tug/Core/Matrix.hpp"
#include <Eigen/src/Core/Array.h> #include <Eigen/src/Core/Array.h>
#include <Eigen/src/Core/AssignEvaluator.h>
#include <Eigen/src/Core/Matrix.h> #include <Eigen/src/Core/Matrix.h>
#include <algorithm> #include <algorithm>
#include <cstddef> #include <cstddef>
@ -190,8 +191,7 @@ private:
cellIndex < 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i); cellIndex < 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i);
// if we have a constant boundary, we use the value of the boundary // if we have a constant boundary, we use the value of the boundary
if (bcElement.getType() == BC_TYPE_CONSTANT) { if (bcElement.getType() == BC_TYPE_CONSTANT) {
currentFlux = currentFlux = currentVelo * bcElement.getValue();
currentVelo * bcElement.getValue() * (innerTimestep) / dx;
continue; continue;
} }
@ -201,8 +201,7 @@ private:
} }
// otherwise we use the concentration of the inflow/outflow cell, // otherwise we use the concentration of the inflow/outflow cell,
// multiplied by the velocity // multiplied by the velocity
currentFlux = currentVelo * oldConcentrations(row_i, cellIndex) * currentFlux = currentVelo * oldConcentrations(row_i, cellIndex);
innerTimestep / dx;
} }
} }
@ -219,8 +218,7 @@ private:
// if we have a constant boundary, we use the value of the // if we have a constant boundary, we use the value of the
// boundary // boundary
if (bcElement.getType() == BC_TYPE_CONSTANT) { if (bcElement.getType() == BC_TYPE_CONSTANT) {
currentFlux = currentFlux = currentVelo * bcElement.getValue();
currentVelo * bcElement.getValue() * (2 * innerTimestep) / dy;
continue; continue;
} }
@ -230,8 +228,7 @@ private:
} }
// otherwise we use the concentration of the inflow/outflow cell, // otherwise we use the concentration of the inflow/outflow cell,
// multiplied by the velocity // multiplied by the velocity
currentFlux = currentVelo * oldConcentrations(cellIndex, col_i) * currentFlux = currentVelo * oldConcentrations(cellIndex, col_i);
innerTimestep / dy;
} }
} }
@ -241,8 +238,10 @@ private:
const T horizontalFlux = const T horizontalFlux =
fluxX(row_i, col_i) - fluxX(row_i, col_i + 1); 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); const T verticalFlux = fluxY(row_i, col_i) - fluxY(row_i + 1, col_i);
newConcentrations(row_i, col_i) = newConcentrations(row_i, col_i) = oldConcentrations(row_i, col_i) +
oldConcentrations(row_i, col_i) + horizontalFlux + verticalFlux; innerTimestep /
porevolume(row_i, col_i) *
(horizontalFlux + verticalFlux);
} }
} }