diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 8997c0e..33aeea2 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -7,6 +7,7 @@ add_executable(BTCS_1D_proto_example BTCS_1D_proto_example.cpp) add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp) add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp) add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp) +add_executable(meeting_example meeting_example.cpp) target_link_libraries(first_example tug) target_link_libraries(second_example tug) target_link_libraries(boundary_example1D tug) @@ -16,6 +17,7 @@ target_link_libraries(BTCS_1D_proto_example tug) target_link_libraries(BTCS_2D_proto_example tug) target_link_libraries(FTCS_2D_proto_example_mdl tug) target_link_libraries(reference-FTCS_2D_closed tug) +target_link_libraries(meeting_example tug) # target_link_libraries(FTCS_2D_proto_example easy_profiler) add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp) diff --git a/examples/meeting_example.cpp b/examples/meeting_example.cpp new file mode 100644 index 0000000..27d5540 --- /dev/null +++ b/examples/meeting_example.cpp @@ -0,0 +1,38 @@ +#include "tug/Boundary.hpp" +#include + +int main(int argc, char *argv[]) { + int row = 20; + int col = 30; + + // grid + Grid grid = Grid(row,col); + MatrixXd concentrations = MatrixXd::Constant(row,col,0); + concentrations(10,4) = 2000; + grid.setConcentrations(concentrations); + MatrixXd alphaX = MatrixXd::Constant(row,col,1); + for (int i = 0; i < row; i++) { + alphaX(i,0) = 0.05; + } + MatrixXd alphaY = MatrixXd::Constant(row,col,1); + grid.setAlpha(alphaX, alphaY); + + // boundary + Boundary bc = Boundary(grid); + bc.setBoundarySideConstant(BC_SIDE_LEFT, 0); + bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0); + bc.setBoundarySideConstant(BC_SIDE_TOP, 0); + bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0); + bc.setBoundaryElementClosed(BC_SIDE_LEFT, 3); + bc.setBoundaryElementClosed(BC_SIDE_LEFT, 4); + bc.setBoundaryElementClosed(BC_SIDE_LEFT, 5); + bc.setBoundaryElementConstant(BC_SIDE_BOTTOM, 17, 100); + + // simulation + Simulation sim = Simulation(grid, bc, BTCS_APPROACH); + sim.setTimestep(0.1); + sim.setIterations(100); + sim.setOutputCSV(CSV_OUTPUT_XTREME); + + sim.run(); +} \ No newline at end of file diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp index 8b7ed21..5e96d94 100644 --- a/src/BTCSv2.cpp +++ b/src/BTCSv2.cpp @@ -342,7 +342,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) { for (int i = 0; i < colMax; i++) { // swap alphas, boundary conditions and sx/sy for column-wise calculation - A = createCoeffMatrix(alphaY, bcLeft, bcRight, rowMax, i, sy); + A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy); b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, bcLeft, bcRight, rowMax, i, sy, sx); row_t1 = solve(A, b); diff --git a/src/Boundary.cpp b/src/Boundary.cpp index 8f248cf..bba7d9d 100644 --- a/src/Boundary.cpp +++ b/src/Boundary.cpp @@ -64,22 +64,36 @@ void Boundary::setBoundarySideClosed(BC_SIDE side) { if(grid.getDim() == 1){ if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ throw_invalid_argument( - "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "For the one-dimensional case, only the BC_SIDE_LEFT and " "BC_SIDE_RIGHT borders exist."); } } - this->boundaries[side] = vector(grid.getRow(), BoundaryElement()); + + int n; + if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { + n = grid.getRow(); + } else { + n = grid.getCol(); + } + this->boundaries[side] = vector(n, BoundaryElement()); } void Boundary::setBoundarySideConstant(BC_SIDE side, double value) { if(grid.getDim() == 1){ if((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)){ throw_invalid_argument( - "For the one-dimensional trap, only the BC_SIDE_LEFT and " + "For the one-dimensional case, only the BC_SIDE_LEFT and " "BC_SIDE_RIGHT borders exist."); } } - this->boundaries[side] = vector(grid.getRow(), BoundaryElement(value)); + + int n; + if (side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT) { + n = grid.getRow(); + } else { + n = grid.getCol(); + } + this->boundaries[side] = vector(n, BoundaryElement(value)); } void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) {