added reference example FTCS 2D closed

examples/CMakeLists.txt

@@ -3,10 +3,12 @@ add_executable(second_example second_example.cpp)
 add_executable(boundary_example1D boundary_example1D.cpp)
 add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.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_2D_proto_example 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) 

examples/reference-FTCS_2D_closed.cpp

@@ -0,0 +1,50 @@
+#include <tug/Simulation.hpp>
+
+int main(int argc, char *argv[]) {
+    int row = 11;
+    int col = 11;
+    int domain_row = 10;
+    int domain_col = 10;
+
+
+    // Grid
+    Grid grid = Grid(row, col);
+    grid.setDomain(domain_row, domain_col);
+
+    MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
+    concentrations(5,5) = 1;
+    grid.setConcentrations(concentrations);
+
+    MatrixXd alpha = MatrixXd::Constant(row, col, 1);
+    for (int i = 0; i < 5; i++) {
+        for (int j = 0; j < 6; j++) {
+            alpha(i, j) = 0.01;
+        }
+    }
+    for (int i = 0; i < 5; i++) {
+        for (int j = 6; j < 11; j++) {
+            alpha(i, j) = 0.001;
+        }
+    }
+    for (int i = 5; i < 11; i++) {
+        for (int j = 6; j < 11; j++) {
+            alpha(i, j) = 0.1;
+        }
+    }
+    grid.setAlpha(alpha, alpha);
+
+
+    // Boundary
+    Boundary bc = Boundary(grid);
+
+
+    // Simulation
+    Simulation sim = Simulation(grid, bc, FTCS_APPROACH);
+    sim.setTimestep(0.001);
+    sim.setIterations(7000);
+    sim.setOutputCSV(CSV_OUTPUT_VERBOSE);
+
+
+    // RUN
+    sim.run();
+}

src/FTCS.cpp

@@ -6,7 +6,7 @@
 using namespace std;
 
 
-double calcAlphaIntercell(double alpha1, double alpha2, bool useHarmonic = false) {
+double calcAlphaIntercell(double alpha1, double alpha2, bool useHarmonic = true) {
     if (useHarmonic) {
         return 2 / ((1/alpha1) + (1/alpha2));
     } else {