diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp index 8cecea9..ce1d5e8 100644 --- a/include/tug/Core/Numeric/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -34,214 +34,6 @@ constexpr T calcChangeInner(T conc_c, T conc_left, T conc_right, T alpha_c, alpha_center_right * conc_right; } -// calculates horizontal change on one cell independent of boundary type -// template -// static inline T calcHorizontalChange(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaX, int &row, -// int &col) { -// return calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) * -// concentrations(row, col + 1) - -// (calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) + -// calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col))) * -// concentrations(row, col) + -// calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col)) * -// concentrations(row, col - 1); -// } -// -// // calculates vertical change on one cell independent of boundary type -// template -// static inline T calcVerticalChange(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaY, int &row, -// int &col) { -// return calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) * -// concentrations(row + 1, col) - -// (calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) + -// calcAlphaIntercell(alphaY(row - 1, col), alphaY(row, col))) * -// concentrations(row, col) + -// calcAlphaIntercell(alphaY(row - 1, col), alphaY(row, col)) * -// concentrations(row - 1, col); -// } - -// calculates horizontal change on one cell with a constant left boundary -// template -// static inline T calcHorizontalChangeLeftBoundaryConstant( -// const RowMajMatMap &concentrations, const RowMajMatMap &alphaX, -// const Boundary &bc, int &row, int &col) { -// return calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) * -// concentrations(row, col + 1) - -// (calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) + -// 2 * alphaX(row, col)) * -// concentrations(row, col) + -// 2 * alphaX(row, col) * bc.getBoundaryElementValue(BC_SIDE_LEFT, -// row); -// } -// -// // calculates horizontal change on one cell with a closed left boundary -// template -// static inline T -// calcHorizontalChangeLeftBoundaryClosed(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaX, int -// &row, int &col) { -// return calcAlphaIntercell(alphaX(row, col + 1), alphaX(row, col)) * -// (concentrations(row, col + 1) - concentrations(row, col)); -// } -// -// // checks boundary condition type for a cell on the left edge of grid -// template -// static inline T -// calcHorizontalChangeLeftBoundary(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaX, Boundary -// &bc, int &row, int &col) { -// if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CONSTANT) { -// return calcHorizontalChangeLeftBoundaryConstant(concentrations, alphaX, -// bc, -// row, col); -// } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CLOSED) -// { -// return calcHorizontalChangeLeftBoundaryClosed(concentrations, alphaX, -// row, -// col); -// } else { -// throw_invalid_argument("Undefined Boundary Condition Type!"); -// } -// } -// -// // calculates horizontal change on one cell with a constant right boundary -// template -// static inline T -// calcHorizontalChangeRightBoundaryConstant(const RowMajMatMap -// &concentrations, -// const RowMajMatMap &alphaX, -// Boundary &bc, int &row, int -// &col) { -// return 2 * alphaX(row, col) * bc.getBoundaryElementValue(BC_SIDE_RIGHT, -// row) - -// (calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col)) + -// 2 * alphaX(row, col)) * -// concentrations(row, col) + -// calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col)) * -// concentrations(row, col - 1); -// } -// -// // calculates horizontal change on one cell with a closed right boundary -// template -// static inline T -// calcHorizontalChangeRightBoundaryClosed(const RowMajMatMap -// &concentrations, -// const RowMajMatMap &alphaX, int -// &row, int &col) { -// return -(calcAlphaIntercell(alphaX(row, col - 1), alphaX(row, col)) * -// (concentrations(row, col) - concentrations(row, col - 1))); -// } -// -// // checks boundary condition type for a cell on the right edge of grid -// template -// static inline T -// calcHorizontalChangeRightBoundary(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaX, -// Boundary &bc, int &row, int &col) { -// if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CONSTANT) { -// return calcHorizontalChangeRightBoundaryConstant(concentrations, alphaX, -// bc, -// row, col); -// } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CLOSED) -// { -// return calcHorizontalChangeRightBoundaryClosed(concentrations, alphaX, -// row, -// col); -// } else { -// throw_invalid_argument("Undefined Boundary Condition Type!"); -// } -// } -// -// // calculates vertical change on one cell with a constant top boundary -// template -// static inline T -// calcVerticalChangeTopBoundaryConstant(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaY, -// Boundary &bc, int &row, int &col) { -// return calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) * -// concentrations(row + 1, col) - -// (calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) + -// 2 * alphaY(row, col)) * -// concentrations(row, col) + -// 2 * alphaY(row, col) * bc.getBoundaryElementValue(BC_SIDE_TOP, col); -// } -// -// // calculates vertical change on one cell with a closed top boundary -// template -// static inline T -// calcVerticalChangeTopBoundaryClosed(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaY, int &row, -// int &col) { -// return calcAlphaIntercell(alphaY(row + 1, col), alphaY(row, col)) * -// (concentrations(row + 1, col) - concentrations(row, col)); -// } -// -// // checks boundary condition type for a cell on the top edge of grid -// template -// static inline T -// calcVerticalChangeTopBoundary(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaY, Boundary &bc, -// int &row, int &col) { -// if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) { -// return calcVerticalChangeTopBoundaryConstant(concentrations, alphaY, bc, -// row, col); -// } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) { -// return calcVerticalChangeTopBoundaryClosed(concentrations, alphaY, row, -// col); -// } else { -// throw_invalid_argument("Undefined Boundary Condition Type!"); -// } -// } -// -// // calculates vertical change on one cell with a constant bottom boundary -// template -// static inline T -// calcVerticalChangeBottomBoundaryConstant(const RowMajMatMap -// &concentrations, -// const RowMajMatMap &alphaY, -// Boundary &bc, int &row, int &col) -// { -// return 2 * alphaY(row, col) * -// bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) - -// (calcAlphaIntercell(alphaY(row, col), alphaY(row - 1, col)) + -// 2 * alphaY(row, col)) * -// concentrations(row, col) + -// calcAlphaIntercell(alphaY(row, col), alphaY(row - 1, col)) * -// concentrations(row - 1, col); -// } -// -// // calculates vertical change on one cell with a closed bottom boundary -// template -// static inline T -// calcVerticalChangeBottomBoundaryClosed(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaY, int -// &row, int &col) { -// return -(calcAlphaIntercell(alphaY(row, col), alphaY(row - 1, col)) * -// (concentrations(row, col) - concentrations(row - 1, col))); -// } -// -// // checks boundary condition type for a cell on the bottom edge of grid -// template -// static inline T -// calcVerticalChangeBottomBoundary(const RowMajMatMap &concentrations, -// const RowMajMatMap &alphaY, Boundary -// &bc, int &row, int &col) { -// if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) { -// return calcVerticalChangeBottomBoundaryConstant(concentrations, alphaY, -// bc, -// row, col); -// } else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == -// BC_TYPE_CLOSED) { -// return calcVerticalChangeBottomBoundaryClosed(concentrations, alphaY, -// row, -// col); -// } else { -// throw_invalid_argument("Undefined Boundary Condition Type!"); -// } -// } - template constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center, T alpha_neighbor, const BoundaryElement &bc) { @@ -428,168 +220,8 @@ static void FTCS_2D(SimulationInput &input, int numThreads) { } } - // inner cells - // these are independent of the boundary condition type - // omp_set_num_threads(10); - // #pragma omp parallel for num_threads(numThreads) - // for (int row = 1; row < rowMax - 1; row++) { - // for (int col = 1; col < colMax - 1; col++) { - // const T &conc_c = concentrations_grid(row, col); - // const T &conc_left = concentrations_grid(row, col - 1); - // const T &conc_right = concentrations_grid(row, col + 1); - // const T &conc_top = concentrations_grid(row + 1, col); - // const T &conc_bottom = concentrations_grid(row - 1, col); - // - // const T &alpha_c = alphaX(row, col); - // const T &alpha_left = alphaX(row, col - 1); - // const T &alpha_right = alphaX(row, col + 1); - // const T &alpha_top = alphaY(row + 1, col); - // const T &alpha_bottom = alphaY(row - 1, col); - // - // const T horizontal_change = calcChangesInner( - // conc_c, conc_left, conc_right, alpha_c, alpha_left, alpha_right); - // - // const T vertical_change = calcChangesInner( - // conc_c, conc_bottom, conc_top, alpha_c, alpha_bottom, alpha_top); - // - // concentrations_t1(row, col) = - // concentrations_grid(row, col) + - // timestep / (deltaRow * deltaRow) * vertical_change + - // timestep / (deltaCol * deltaCol) * horizontal_change; - // } - // } - - // boundary conditions - // left without corners / looping over rows - // hold column constant at index 0 - // int col = 0; - // #pragma omp parallel for num_threads(numThreads) - // for (int row = 1; row < rowMax - 1; row++) { - // const T horizontal_change = calcChangeBoundary( - // concentrations_grid(row, col), concentrations_grid(row, col + 1), - // alphaX(row, col), alphaX(row, col + 1), - // bc.getBoundaryElement(BC_SIDE_LEFT, row)); - // - // const T vertical_change = calcChangesInner( - // concentrations_grid(row, col), concentrations_grid(row - 1, col), - // concentrations_grid(row + 1, col), alphaX(row, col), - // alphaY(row - 1, col), alphaY(row + 1, col)); - // - // concentrations_t1(row, col) = - // concentrations_grid(row, col) + - // timestep / (deltaCol * deltaCol) * horizontal_change + - // timestep / (deltaRow * deltaRow) * vertical_change; - // } - // - // // right without corners / looping over rows - // // hold column constant at max index - // col = colMax - 1; - // #pragma omp parallel for num_threads(numThreads) - // for (int row = 1; row < rowMax - 1; row++) { - // - // concentrations_t1(row, col) = - // concentrations_grid(row, col) + - // timestep / (deltaCol * deltaCol) * - // (calcHorizontalChangeRightBoundary(concentrations_grid, alphaX, - // bc, - // row, col)) + - // timestep / (deltaRow * deltaRow) * - // (calcVerticalChange(concentrations_grid, alphaY, row, col)); - // } - // - // // top without corners / looping over columns - // // hold row constant at index 0 - // int row = 0; - // #pragma omp parallel for num_threads(numThreads) - // for (int col = 1; col < colMax - 1; col++) { - // concentrations_t1(row, col) = - // concentrations_grid(row, col) + - // timestep / (deltaRow * deltaRow) * - // (calcVerticalChangeTopBoundary(concentrations_grid, alphaY, bc, - // row, - // col)) + - // timestep / (deltaCol * deltaCol) * - // (calcHorizontalChange(concentrations_grid, alphaX, row, col)); - // } - // - // // bottom without corners / looping over columns - // // hold row constant at max index - // row = rowMax - 1; - // #pragma omp parallel for num_threads(numThreads) - // for (int col = 1; col < colMax - 1; col++) { - // concentrations_t1(row, col) = - // concentrations_grid(row, col) + - // timestep / (deltaRow * deltaRow) * - // (calcVerticalChangeBottomBoundary(concentrations_grid, alphaY, - // bc, - // row, col)) + - // timestep / (deltaCol * deltaCol) * - // (calcHorizontalChange(concentrations_grid, alphaX, row, col)); - // } - // - // // corner top left - // // hold row and column constant at 0 - // row = 0; - // col = 0; - // concentrations_t1(row, col) = - // concentrations_grid(row, col) + - // timestep / (deltaCol * deltaCol) * - // (calcHorizontalChangeLeftBoundary(concentrations_grid, alphaX, - // bc, - // row, col)) + - // timestep / (deltaRow * deltaRow) * - // (calcVerticalChangeTopBoundary(concentrations_grid, alphaY, bc, - // row, - // col)); - // - // // corner top right - // // hold row constant at 0 and column constant at max index - // row = 0; - // col = colMax - 1; - // concentrations_t1(row, col) = - // concentrations_grid(row, col) + - // timestep / (deltaCol * deltaCol) * - // (calcHorizontalChangeRightBoundary(concentrations_grid, alphaX, - // bc, - // row, col)) + - // timestep / (deltaRow * deltaRow) * - // (calcVerticalChangeTopBoundary(concentrations_grid, alphaY, bc, - // row, - // col)); - // - // // corner bottom left - // // hold row constant at max index and column constant at 0 - // row = rowMax - 1; - // col = 0; - // concentrations_t1(row, col) = - // concentrations_grid(row, col) + - // timestep / (deltaCol * deltaCol) * - // (calcHorizontalChangeLeftBoundary(concentrations_grid, alphaX, - // bc, - // row, col)) + - // timestep / (deltaRow * deltaRow) * - // (calcVerticalChangeBottomBoundary(concentrations_grid, alphaY, - // bc, - // row, col)); - // - // // corner bottom right - // // hold row and column constant at max index - // row = rowMax - 1; - // col = colMax - 1; - // concentrations_t1(row, col) = - // concentrations_grid(row, col) + - // timestep / (deltaCol * deltaCol) * - // (calcHorizontalChangeRightBoundary(concentrations_grid, alphaX, - // bc, - // row, col)) + - // timestep / (deltaRow * deltaRow) * - // (calcVerticalChangeBottomBoundary(concentrations_grid, alphaY, - // bc, - // row, col)); - // overwrite obsolete concentrations concentrations_grid = concentrations_t1; - // } } // entry point; differentiate between 1D and 2D grid