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

refactor: implement coeff boundary functions as template constexpr

Browse Source
This commit is contained in:
Max Lübke 2023-09-14 16:21:45 +02:00
parent ef1ccd4c14
commit 6c1ccb90fd
2 changed files with 24 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

69
src/BTCS.cpp
View File

@ -14,6 +14,7 @@
#include <cstddef> #include <cstddef>
#include <tug/Boundary.hpp> #include <tug/Boundary.hpp>
#include <tug/Grid.hpp> #include <tug/Grid.hpp>
#include <utility>
#ifdef _OPENMP #ifdef _OPENMP
#include <omp.h> #include <omp.h>
@ -21,52 +22,26 @@
#define omp_get_thread_num() 0 #define omp_get_thread_num() 0
#endif #endif
// calculates coefficient for left boundary in constant case // calculates coefficient for boundary in constant case
static inline std::tuple<double, double> template <class T>
calcLeftBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, double sx) { constexpr std::pair<T, T> calcBoundaryCoeffConstant(T alpha_center,
const double centerCoeff = T alpha_side, T sx) {
1 + sx * (calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)) + const T centerCoeff = 1 + sx * (calcAlphaIntercell(alpha_center, alpha_side) +
2 * alpha(rowIndex, 0)); 2 * alpha_center);
const double rightCoeff = const T sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side);
-sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1));
return {centerCoeff, rightCoeff}; return {centerCoeff, sideCoeff};
} }
// calculates coefficient for left boundary in closed case // calculates coefficient for boundary in closed case
static inline std::tuple<double, double> template <class T>
calcLeftBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, double sx) { constexpr std::pair<T, T> calcBoundaryCoeffClosed(T alpha_center, T alpha_side,
const double centerCoeff = T sx) {
1 + sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1)); const T centerCoeff = 1 + sx * calcAlphaIntercell(alpha_center, alpha_side);
const double rightCoeff =
-sx * calcAlphaIntercell(alpha(rowIndex, 0), alpha(rowIndex, 1));
return {centerCoeff, rightCoeff}; const T sideCoeff = -sx * calcAlphaIntercell(alpha_center, alpha_side);
}
// calculates coefficient for right boundary in constant case return {centerCoeff, sideCoeff};
static inline std::tuple<double, double>
calcRightBoundaryCoeffConstant(Eigen::MatrixXd &alpha, int rowIndex, int n,
double sx) {
const double leftCoeff =
-sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
const double centerCoeff =
1 + sx * (calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n)) +
2 * alpha(rowIndex, n));
return {leftCoeff, centerCoeff};
}
// calculates coefficient for right boundary in closed case
static inline std::tuple<double, double>
calcRightBoundaryCoeffClosed(Eigen::MatrixXd &alpha, int rowIndex, int n,
double sx) {
const double leftCoeff =
-sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
const double centerCoeff =
1 + sx * calcAlphaIntercell(alpha(rowIndex, n - 1), alpha(rowIndex, n));
return {leftCoeff, centerCoeff};
} }
// creates coefficient matrix for next time step from alphas in x-direction // creates coefficient matrix for next time step from alphas in x-direction
@ -84,14 +59,14 @@ createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector<BoundaryElement> &bcLeft,
switch (bcLeft[rowIndex].getType()) { switch (bcLeft[rowIndex].getType()) {
case BC_TYPE_CONSTANT: { case BC_TYPE_CONSTANT: {
auto [centerCoeffTop, rightCoeffTop] = auto [centerCoeffTop, rightCoeffTop] =
calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx); calcBoundaryCoeffConstant(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
cm.insert(0, 0) = centerCoeffTop; cm.insert(0, 0) = centerCoeffTop;
cm.insert(0, 1) = rightCoeffTop; cm.insert(0, 1) = rightCoeffTop;
break; break;
} }
case BC_TYPE_CLOSED: { case BC_TYPE_CLOSED: {
auto [centerCoeffTop, rightCoeffTop] = auto [centerCoeffTop, rightCoeffTop] =
calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx); calcBoundaryCoeffClosed(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
cm.insert(0, 0) = centerCoeffTop; cm.insert(0, 0) = centerCoeffTop;
cm.insert(0, 1) = rightCoeffTop; cm.insert(0, 1) = rightCoeffTop;
break; break;
@ -119,15 +94,15 @@ createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector<BoundaryElement> &bcLeft,
switch (bcRight[rowIndex].getType()) { switch (bcRight[rowIndex].getType()) {
case BC_TYPE_CONSTANT: { case BC_TYPE_CONSTANT: {
auto [leftCoeffBottom, centerCoeffBottom] = auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffConstant(
calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx); alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
cm.insert(n, n - 1) = leftCoeffBottom; cm.insert(n, n - 1) = leftCoeffBottom;
cm.insert(n, n) = centerCoeffBottom; cm.insert(n, n) = centerCoeffBottom;
break; break;
} }
case BC_TYPE_CLOSED: { case BC_TYPE_CLOSED: {
auto [leftCoeffBottom, centerCoeffBottom] = auto [centerCoeffBottom, leftCoeffBottom] =
calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx); calcBoundaryCoeffClosed(alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
cm.insert(n, n - 1) = leftCoeffBottom; cm.insert(n, n - 1) = leftCoeffBottom;
cm.insert(n, n) = centerCoeffBottom; cm.insert(n, n) = centerCoeffBottom;
break; break;

3
src/TugUtils.hpp
View File

@ -25,7 +25,8 @@
}) })
// calculates arithmetic or harmonic mean of alpha between two cells // calculates arithmetic or harmonic mean of alpha between two cells
constexpr double calcAlphaIntercell(double alpha1, double alpha2, template <typename T>
constexpr T calcAlphaIntercell(T alpha1, T alpha2,
bool useHarmonic = true) { bool useHarmonic = true) {
if (useHarmonic) { if (useHarmonic) {
return double(2) / ((double(1) / alpha1) + (double(1) / alpha2)); return double(2) / ((double(1) / alpha1) + (double(1) / alpha2));