implemented some util methods and started with a first test case

examples/reference-FTCS_2D_closed.cpp

@@ -1,5 +1,4 @@
 #include <tug/Simulation.hpp>
-#include "csv2eigen.cpp"
 #include "Eigen/Core"
 
 int main(int argc, char *argv[]) {
@@ -42,17 +41,13 @@ int main(int argc, char *argv[]) {
 
     // Simulation
     Simulation sim = Simulation(grid, bc, FTCS_APPROACH);
-    //sim.setTimestep(0.001);
-    sim.setIterations(2);
+    sim.setTimestep(0.001);
+    sim.setIterations(7000);
     sim.setOutputCSV(CSV_OUTPUT_ON);
+    sim.setOutputConsole(CONSOLE_OUTPUT_ON);
 
-    MatrixXd mymatrix = CSV2Eigen("/Users/hannessigner/Documents/tug/build/examples/FTCS_11_11_2-1.csv");
-
-    cout << "Matrix start:" << endl;
-    cout << mymatrix << endl;
-
-    //bool r = grid.isApprox(mymatrix);
 
     // RUN
     sim.run();
+
 }

src/Simulation.cpp

@@ -25,12 +25,13 @@ Simulation::Simulation(Grid grid, Boundary bc, APPROACH approach) : grid(grid),
     double maxAlphaY = grid.getAlphaY().maxCoeff();
     double maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY;
 
-    //double maxStableTimestep = minDelta / (2*maxAlpha); // Formula from Marco --> seems to be unstable
+    double maxStableTimestepMdl = minDelta / (2*maxAlpha); // Formula from Marco --> seems to be unstable
     double maxStableTimestep = 1 / (4 * maxAlpha * ((1/deltaRowSquare) + (1/deltaColSquare))); // Formula from Wikipedia
 
-    cout << maxStableTimestep << endl;
+    // cout << "Max stable time step MDL: " << maxStableTimestepMdl << endl;
+    // cout << "Max stable time step: " << maxStableTimestep << endl;
 
-    this->timestep = maxStableTimestep;
+    this->timestep = maxStableTimestep; 
 
     
     this->iterations = 1000;

test/CMakeLists.txt

@@ -11,7 +11,7 @@ if(NOT DOCTEST_LIB)
   FetchContent_MakeAvailable(DocTest)
 endif()
 
-add_executable(testTug setup.cpp testBoundaryCondition.cpp testDiffusion.cpp)
+add_executable(testTug setup.cpp testBoundaryCondition.cpp testDiffusion.cpp testSimulation.cpp)
 target_link_libraries(testTug doctest tug)
 
 add_custom_target(

test/FTCS_11_11_7000.csv

@@ -0,0 +1,13 @@
+0 0 0 0 0 0 0 0 0 0 0
+1.88664e-08 3.39962e-08 7.57021e-08 1.76412e-07 4.15752e-07 9.00973e-07 3.65403e-09 9.6579e-12 3.59442e-13 3.42591e-14 3.27595e-15
+1.19102e-06 1.95195e-06 3.92165e-06 8.30575e-06 1.78976e-05 3.60742e-05 2.02843e-07 2.24659e-09 2.35085e-10 2.64988e-11 2.90933e-12
+5.85009e-05 8.57948e-05 0.000151499 0.000284105 0.00054607 0.00100251 1.18494e-05 8.26706e-07 1.26394e-07 1.70309e-08 2.20525e-09
+0.00202345 0.00258511 0.00381783 0.00599829 0.00972689 0.0154873 0.0011152 0.000247309 5.10506e-05 8.64727e-06 1.37747e-06
+0.0205848 0.0217651 0.0238282 0.0262762 0.0285812 0.0303808 0.0374255 0.0204234 0.00674813 0.00160264 0.000338852
+0.0199112 0.0210265 0.0229587 0.0252019 0.0272007 0.0285225 0.0310896 0.0156681 0.00495399 0.00114176 0.000235573
+0.0184589 0.0194561 0.0211596 0.0230704 0.0246215 0.0253266 0.0228374 0.010038 0.00292986 0.000638799 0.000125942
+0.0166888 0.0175517 0.0190015 0.0205611 0.0216793 0.0218694 0.0160278 0.00593307 0.00155164 0.000312583 5.76964e-05
+0.0151262 0.0158758 0.0171155 0.0183949 0.019193 0.0190523 0.0115557 0.00356455 0.000817461 0.000148892 2.51893e-05
+0.0142177 0.0149034 0.0160255 0.0171522 0.0177843 0.0174891 0.0093846 0.0025221 0.000515469 8.51039e-05 1.31328e-05
+
+

test/TestUtils.cpp

@@ -1,15 +1,22 @@
+#include <ios>
 #include <iostream>
 #include <Eigen/Core>
+#include <Eigen/Dense>
 #include <fstream>
+#include <sstream>
+#include <stdexcept>
 
 using namespace std;
 using namespace Eigen;
 
-MatrixXd CSV2Eigen(string file2Convert){
+MatrixXd CSV2Eigen(string file2Convert) {
 
     vector<double> matrixEntries;
 
     ifstream matrixDataFile(file2Convert);
+    if (matrixDataFile.fail()) {
+        throw invalid_argument("File probably non-existent!");
+    }
 
     string matrixRowString;
     string matrixEntry;
@@ -20,9 +27,21 @@ MatrixXd CSV2Eigen(string file2Convert){
         while(getline(matrixRowStringStream, matrixEntry, ' ')){
             matrixEntries.push_back(stod(matrixEntry));
         }
-        matrixRowNumber++;
+        if (matrixRowString.length() > 1) {
+            matrixRowNumber++;
+        }
     }
 
     return Map<Matrix<double, Dynamic, Dynamic, RowMajor>>(matrixEntries.data(), matrixRowNumber, matrixEntries.size() / matrixRowNumber);
+}
 
+bool checkSimilarity(MatrixXd a, MatrixXd b, double precision=1e-5) {
+    return a.isApprox(b, precision);
+}
+
+bool checkSimilarityV2(MatrixXd a, MatrixXd b, double maxDiff) {
+
+    MatrixXd diff = a - b;
+    double maxCoeff = diff.maxCoeff();
+    return maxCoeff < maxDiff;
 }

test/testSimulation.cpp

@@ -0,0 +1,66 @@
+#include <doctest/doctest.h>
+#include <tug/Simulation.hpp>
+
+#include "TestUtils.cpp"
+
+static Grid setupSimulation() {
+    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_ON);
+    // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
+
+
+    // RUN
+    sim.run();
+
+    return grid;
+}
+
+TEST_CASE("equality to reference matrix") {
+    MatrixXd reference = CSV2Eigen("/Users/philipp/forschungsprojekt/tug/test/FTCS_11_11_7000.csv");
+
+    Grid grid = setupSimulation();
+
+    cout << reference << endl << endl;
+    cout << grid.getConcentrations() << endl;
+    CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true);
+}