working BTCS 1D and 2D simulation with constant boundary

examples/BTCS_1D_proto_example.cpp

@@ -0,0 +1,48 @@
+#include <tug/Simulation.hpp>
+
+int main(int argc, char *argv[]) {
+    // **************
+    // **** GRID ****
+    // **************
+
+    // create a linear grid with 20 cells
+    int cells = 20;
+    Grid grid = Grid(cells);
+
+    MatrixXd concentrations = MatrixXd::Constant(1,20,0);
+    concentrations(0,0) = 2000;
+    // TODO add option to set concentrations with a vector in 1D case
+    grid.setConcentrations(concentrations);
+
+
+    // ******************
+    // **** BOUNDARY ****
+    // ******************
+
+    // create a boundary with constant values
+    Boundary bc = Boundary(grid);
+    bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
+    bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+
+
+    // ************************
+    // **** SIMULATION ENV ****
+    // ************************
+
+    // set up a simulation environment
+    Simulation simulation = Simulation(grid, bc, BTCS_APPROACH); // grid,boundary,simulation-approach
+
+    // set the timestep of the simulation
+    simulation.setTimestep(0.1); // timestep
+
+    // set the number of iterations
+    simulation.setIterations(100);
+
+    // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
+    simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
+    
+    // **** RUN SIMULATION ****
+
+    // run the simulation
+    simulation.run();
+}

examples/CMakeLists.txt

@@ -1,18 +1,20 @@
 add_executable(first_example first_example.cpp)
 add_executable(second_example second_example.cpp)
 add_executable(boundary_example1D boundary_example1D.cpp)
+add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp)
 add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp)
+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(FTCS_1D_proto_example FTCS_1D_proto_example.cpp)
 add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp)
 target_link_libraries(first_example tug)
 target_link_libraries(second_example tug)
 target_link_libraries(boundary_example1D tug)
+target_link_libraries(FTCS_1D_proto_example tug)
 target_link_libraries(FTCS_2D_proto_example tug)
+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(FTCS_1D_proto_example tug)
 target_link_libraries(reference-FTCS_2D_closed tug)
 # target_link_libraries(FTCS_2D_proto_example easy_profiler) 
 

src/BTCSv2.cpp

@@ -6,11 +6,11 @@
 using namespace Eigen;
 
 
-static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, int dim, int rowIndex, double sx) {
+static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, int numCols, int rowIndex, double sx) {
 
     // square matrix of column^2 dimension for the coefficients
-    SparseMatrix<double> cm(dim, dim);
-    cm.reserve(VectorXi::Constant(dim, 3));
+    SparseMatrix<double> cm(numCols, numCols);
+    cm.reserve(VectorXi::Constant(numCols, 3));
 
     // left column
     cm.insert(0,0) = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1))
@@ -18,13 +18,14 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, int dim, int rowI
     cm.insert(0,1) = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
 
     // inner columns
-    int n = dim-1;
+    int n = numCols-1;
     for (int i = 1; i < n; i++) {
         cm.insert(i,i-1) = -sx * calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i));
         cm.insert(i,i) = 1 + sx * (
                             calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1))
-                            + calcAlphaIntercell(alpha(rowIndex,i-1), alpha(i,1))
-                            );
+                            + calcAlphaIntercell(alpha(rowIndex,i-1), alpha(rowIndex,i))
+                            )
+                        ;
         cm.insert(i,i+1) = -sx * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex,i+1));
     }
 
@@ -96,21 +97,15 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
     }
 
     // first column -> additional fixed concentration change from perpendicular dimension
-    sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft(0);
+    sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft(rowIndex);
 
     // last column -> additional fixed concentration change from perpendicular dimension
-    sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight(length-1);
+    sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight(rowIndex);
 
     return sv;
 }
 
 
-// BTCS solution for 1D grid
-static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
-    // TODO
-}
-
-
 static VectorXd solve(SparseMatrix<double> &A, VectorXd &b) {
     SparseLU<SparseMatrix<double>> solver;
     solver.analyzePattern(A);
@@ -120,6 +115,38 @@ static VectorXd solve(SparseMatrix<double> &A, VectorXd &b) {
 }
 
 
+// BTCS solution for 1D grid 
+static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
+    int length = grid.getLength();
+    double sx = timestep / (grid.getDeltaCol() * grid.getDeltaCol()); // TODO create method getDelta() for 1D case
+
+    VectorXd concentrations_t1(length);
+
+    SparseMatrix<double> A;
+    VectorXd b(length);
+
+    MatrixXd alpha = grid.getAlpha();
+    VectorXd bcLeft = bc.getBoundarySideValues(BC_SIDE_LEFT);
+    VectorXd bcRight = bc.getBoundarySideValues(BC_SIDE_RIGHT);
+
+    MatrixXd concentrations = grid.getConcentrations();
+    A = createCoeffMatrix(alpha, length, 0, sx); // this is exactly same as in 2D
+    for (int i = 0; i < length; i++) {
+        b(i) = concentrations(0,i); // TODO check if this is really the only thing on right hand side of equation in 1D?
+    }
+    b(0) += 2 * sx * alpha(0,0) * bcLeft(0);
+    b(length-1) += 2 * sx * alpha(0,length-1) * bcRight(0);
+
+    concentrations_t1 = solve(A, b);
+
+    for (int j = 0; j < length; j++) {
+        concentrations(0,j) = concentrations_t1(j);
+    }
+
+    grid.setConcentrations(concentrations);
+}
+
+
 // BTCS solution for 2D grid
 static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
     int rowMax = grid.getRow();