Merge branch 'hannes-philipp' of git.gfz-potsdam.de:naaice/tug into hannes-philipp

2025-12-13 09:28:23 +01:00 · 2023-08-29 10:42:23 +02:00 · 2023-08-29 10:42:23 +02:00 · 8fcc8812e7
commit 8fcc8812e7
parent a0ab764870 4cb51f4241
15 changed files with 684 additions and 126 deletions
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@ -4,6 +4,7 @@ stages:
    - build
    - test
    - static_analyze
+    - doc

 build_release:
  stage: build
@ -22,6 +23,23 @@ test:
    script:
        - ./build/test/testTug

+pages:
+  stage: doc
+  image: python:slim
+  before_script:
+    - apt-get update && apt-get install --no-install-recommends -y graphviz imagemagick doxygen make
+    - pip install --upgrade pip && pip install Sphinx Pillow breathe sphinx-rtd-theme
+    - mkdir public
+  script:
+    - pushd docs_sphinx
+    - make html
+    - popd && mv docs_sphinx/_build/html/* public/
+  artifacts:
+    paths:
+      - public
+  rules:
+    - if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH
+
 lint:
  stage: static_analyze
  before_script:
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@ -1,32 +0,0 @@
-# .readthedocs.yaml
-# Read the Docs configuration file
-# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
-
-# Required
-version: 2
-
-# Set the OS, Python version and other tools you might need
-build:
-  os: ubuntu-22.04
-  tools:
-    python: "3.11"
-    # You can also specify other tool versions:
-    # nodejs: "19"
-    # rust: "1.64"
-    # golang: "1.19"
-
-# Build documentation in the "docs/" directory with Sphinx
-sphinx:
-   configuration: docs_sphinx/conf.py
-
-# Optionally build your docs in additional formats such as PDF and ePub
-formats:
-   - pdf
-#    - epub
-
-# Optional but recommended, declare the Python requirements required
-# to build your documentation
-# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html
-python:
-   install:
-   - requirements: docs_sphinx/requirements.txt
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -7,7 +7,8 @@ set(CMAKE_CXX_STANDARD 17)

 find_package(Eigen3 REQUIRED NO_MODULE)
 find_package(OpenMP)
-find_package(easy_profiler)
+# find_package(easy_profiler)
+# option(EASY_OPTION_LOG "Verbose easy_profiler" 1)

 ## SET(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -O2 -mfma")
 option(TUG_USE_OPENMP "Compile with OpenMP support" ON)
@ -40,4 +41,4 @@ if(TUG_ENABLE_TESTING)
  add_subdirectory(test)
 endif()

-add_subdirectory(examples)
+add_subdirectory(examples)
--- a/doc/ADI_scheme.org
+++ b/doc/ADI_scheme.org
@ -1,6 +1,7 @@
 #+TITLE: Finite Difference Schemes for the numerical solution of heterogeneous diffusion equation in 2D
 #+LaTeX_CLASS_OPTIONS: [a4paper,10pt]
 #+LATEX_HEADER: \usepackage{fullpage}
+#+LATEX_HEADER: \usepackage{charter}
 #+LATEX_HEADER: \usepackage{amsmath, systeme, cancel, xcolor}
 #+OPTIONS: toc:nil

--- a/doc/ADI_scheme.pdf
+++ b/doc/ADI_scheme.pdf
--- a/doc/ValidationHetDiff.org
+++ b/doc/ValidationHetDiff.org
@ -1,12 +1,12 @@
-#+TITLE: 2D Validation Examples
+#+TITLE: Validation Examples for 2D Heterogeneous Diffusion
 #+AUTHOR: MDL <delucia@gfz-potsdam.de>
-#+DATE: 2023-07-31
+#+DATE: 2023-08-26
 #+STARTUP: inlineimages
 #+LATEX_CLASS_OPTIONS: [a4paper,9pt]
 #+LATEX_HEADER: \usepackage{fullpage}
 #+LATEX_HEADER: \usepackage{amsmath, systeme}
 #+LATEX_HEADER: \usepackage{graphicx}
-#+LATEX_HEADER: \usepackage{}
+#+LATEX_HEADER: \usepackage{charter}
 #+OPTIONS: toc:nil


@ -38,7 +38,7 @@ constant in 4 quadrants:
 The relevant part of the R script used to produce these results is
 presented in listing 1; the whole script is at [[file:scripts/HetDiff.R]].
 A visualization of the output of the reference simulation is given in
-figure [[#fig:1][1]].
+figure [[fig:1][1]].

 Note: all results from this script are stored in the =outc= matrix by
 the =deSolve= function. I stored a different version into
@ -47,68 +47,69 @@ for each time step including initial conditions) and 121 rows, one for
 each domain element, with no headers.

 #+caption: Result of ReacTran/deSolve solution of the above problem at 4
+#+name: fig:1
 [[./images/deSolve_AlphaHet1.png]]


 #+name: lst:1
 #+begin_src R :language R :frame single :caption Listing 1, generate reference simulation using R packages deSolve/ReacTran :captionpos b :label lst:1
-library(ReacTran)
-library(deSolve)
+  library(ReacTran)
+  library(deSolve)

-## harmonic mean
-harm <- function(x,y) {
-  if (length(x) != 1 || length(y) != 1)
-  stop("x & y have different lengths")
-  2/(1/x+1/y)
-}
-
-N     <- 11          # number of grid cells
-ini   <- 1           # initial value at x=0
-N2    <- ceiling(N/2)
-L     <- 10          # domain side
-
-## Define diff.coeff per cell, in  4 quadrants
-alphas <- matrix(0, N, N) 
-alphas[1:N2, 1:N2] <- 1
-alphas[1:N2, seq(N2+1,N)] <- 0.1
-alphas[seq(N2+1,N), 1:N2] <- 0.01
-alphas[seq(N2+1,N), seq(N2+1,N)] <- 0.001
-
-cmpharm <- function(x) {
-  y <- c(0, x, 0)
-  ret <- numeric(length(x)+1)
-  for (i in seq(2, length(y))) {
-    ret[i-1] <- harm(y[i], y[i-1])
+  ## harmonic mean
+  harm <- function(x,y) {
+    if (length(x) != 1 || length(y) != 1)
+    stop("x & y have different lengths")
+    2/(1/x+1/y)
  }
-  ret
-}

-## Construction of the 2D grid
-x.grid <- setup.grid.1D(x.up = 0, L = L, N = N)
-y.grid <- setup.grid.1D(x.up = 0, L = L, N = N)
-grid2D <- setup.grid.2D(x.grid, y.grid)
-dx <- dy <- L/N
+  N     <- 11          # number of grid cells
+  ini   <- 1           # initial value at x=0
+  N2    <- ceiling(N/2)
+  L     <- 10          # domain side

-D.grid <- list()
-## Diffusion coefs on x-interfaces
-D.grid$x.int <- apply(alphas, 1, cmpharm)
-## Diffusion coefs on y-interfaces
-D.grid$y.int <- t(apply(alphas, 2, cmpharm))
+  ## Define diff.coeff per cell, in  4 quadrants
+  alphas <- matrix(0, N, N) 
+  alphas[1:N2, 1:N2] <- 1
+  alphas[1:N2, seq(N2+1,N)] <- 0.1
+  alphas[seq(N2+1,N), 1:N2] <- 0.01
+  alphas[seq(N2+1,N), seq(N2+1,N)] <- 0.001

-# The model
-Diff2Dc <- function(t, y, parms) {
-  CONC  <- matrix(nrow = N, ncol = N, data = y)
-  dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC
-  return(list(dCONC))
-}
+  cmpharm <- function(x) {
+    y <- c(0, x, 0)
+    ret <- numeric(length(x)+1)
+    for (i in seq(2, length(y))) {
+      ret[i-1] <- harm(y[i], y[i-1])
+    }
+    ret
+  }

-## initial condition: 0 everywhere, except in central point
-y <- matrix(nrow = N, ncol = N, data = 0)
-y[N2, N2] <- ini  # initial concentration in the central point...
+  ## Construction of the 2D grid
+  x.grid <- setup.grid.1D(x.up = 0, L = L, N = N)
+  y.grid <- setup.grid.1D(x.up = 0, L = L, N = N)
+  grid2D <- setup.grid.2D(x.grid, y.grid)
+  dx <- dy <- L/N

-## solve for 10 time units
-times <- 0:10
-outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL,
-               dim = c(N, N), lrw = 1860000)
+  D.grid <- list()
+  ## Diffusion coefs on x-interfaces
+  D.grid$x.int <- apply(alphas, 1, cmpharm)
+  ## Diffusion coefs on y-interfaces
+  D.grid$y.int <- t(apply(alphas, 2, cmpharm))
+
+  # The model
+  Diff2Dc <- function(t, y, parms) {
+    CONC  <- matrix(nrow = N, ncol = N, data = y)
+    dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC
+    return(list(dCONC))
+  }
+
+  ## initial condition: 0 everywhere, except in central point
+  y <- matrix(nrow = N, ncol = N, data = 0)
+  y[N2, N2] <- ini  # initial concentration in the central point...
+
+  ## solve for 10 time units
+  times <- 0:10
+  outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL,
+		 dim = c(N, N), lrw = 1860000)
 #+end_src

--- a/docs_sphinx/index.rst
+++ b/docs_sphinx/index.rst
@ -5,12 +5,26 @@

 Welcome to Tug's documentation!
 ===============================
-Welcome to the documentation of the TUG project, a simulation program 
-for solving one- and two-dimensional diffusion problems with heterogeneous diffusion coefficients, more 
-generally, for solving the following differential equation
+
+Welcome to the documentation of the TUG project, a simulation program
+for solving transport equations in one- and two-dimensional uniform
+grids using cell centered finite differences.
+
+---------
+Diffusion
+---------
+
+TUG can solve diffusion problems with heterogeneous and anisotropic
+diffusion coefficients. The partial differential equation expressing
+diffusion reads:

 .. math::
-   \frac{\partial C}{\partial t} = \alpha_x \frac{\partial^2 C}{\partial x^2} + \alpha_y \frac{\partial^2 C}{\partial y^2}.
+   \frac{\partial C}{\partial t} =  \nabla \cdot \left[ \mathbf{\alpha} \nabla C \right]
+
+In 2D, the equation reads:
+   
+.. math::
+   \frac{\partial C}{\partial t} =  \frac{\partial}{\partial x}\left[ \alpha_x \frac{\partial C}{\partial x}\right] + \frac{\partial}{\partial y}\left[ \alpha_y \frac{\partial C}{\partial y}\right]

 .. toctree::
   :maxdepth: 2
--- a/examples/BTCS_2D_proto_example.cpp
+++ b/examples/BTCS_2D_proto_example.cpp
@ -8,8 +8,8 @@ int main(int argc, char *argv[]) {
    // **************
    // profiler::startListen();
    // create a grid with a 20 x 20 field
-    int row = 20;
-    int col = 20;
+    int row = 40;
+    int col = 50;
    Grid grid = Grid(row,col);

    // (optional) set the domain, e.g.:
@ -19,7 +19,7 @@ int main(int argc, char *argv[]) {
    // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value
    // grid.setConcentrations(concentrations);
    MatrixXd concentrations = MatrixXd::Constant(row,col,0);
-    concentrations(0,0) = 2000;
+    concentrations(10,10) = 2000;
    grid.setConcentrations(concentrations);
    

@ -67,7 +67,7 @@ int main(int argc, char *argv[]) {
    simulation.setTimestep(0.1); // timestep

    // set the number of iterations
-    simulation.setIterations(100);
+    simulation.setIterations(300);

    // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
    simulation.setOutputCSV(CSV_OUTPUT_XTREME);
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@ -4,6 +4,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(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp)
 add_executable(profiling_openmp profiling_openmp.cpp)
+
 target_link_libraries(FTCS_1D_proto_example tug)
 target_link_libraries(FTCS_2D_proto_example tug)
 target_link_libraries(BTCS_1D_proto_example tug)
--- a/include/tug/Boundary.hpp
+++ b/include/tug/Boundary.hpp
@ -98,8 +98,6 @@ class BoundaryElement {


 /**
- * This class implements the functionality and management of the boundary
- * conditions in the grid to be simulated.
 * This class implements the functionality and management of the boundary
 * conditions in the grid to be simulated.
 */
--- a/include/tug/Grid.hpp
+++ b/include/tug/Grid.hpp
@ -10,6 +10,7 @@

 using namespace Eigen;

+// TODO document default values and perhaps adjust them
 class Grid {
    public:

--- a/include/tug/Simulation.hpp
+++ b/include/tug/Simulation.hpp
@ -17,6 +17,7 @@ using namespace std;
 enum APPROACH {
    FTCS_APPROACH, // Forward Time-Centered Space
    BTCS_APPROACH, // Backward Time-Centered Space solved with EigenLU solver
+    CRANK_NICOLSON_APPROACH 
 };

 /**
@ -67,8 +68,10 @@ enum TIME_MEASURE {
 class Simulation {
    public:
      /**
-       * @brief Set up a runnable simulation environment with the largest stable
-       *        time step and 1000 iterations by passing the required parameters.
+       * @brief Set up a simulation environment. The timestep and number of iterations
+       *        must be set. For the BTCS approach, the Thomas algorithm is used as 
+       *        the default linear equation solver as this is faster for tridiagonal
+       *        matrices. CSV output, console output and time measure are off by default. 
       *
       * @param grid Valid grid object
       * @param bc Valid boundary condition object
@ -77,7 +80,7 @@ class Simulation {
      Simulation(Grid &grid, Boundary &bc, APPROACH approach);

      /**
-       * @brief Set the option to output the results to a CSV file.
+       * @brief Set the option to output the results to a CSV file. Off by default.
       *
       *
       * @param csv_output Valid output option. The following options can be set
@ -93,7 +96,7 @@ class Simulation {
      void setOutputCSV(CSV_OUTPUT csv_output);

      /**
-       * @brief Set the options for outputting information to the console.
+       * @brief Set the options for outputting information to the console. Off by default. 
       *
       * @param console_output Valid output option. The following options can be set
       *                       here:
@ -103,16 +106,17 @@ class Simulation {
       */
      void setOutputConsole(CONSOLE_OUTPUT console_output);

+    // TODO document method
      /**
-       * @brief Set the Time Measure object
+       * @brief Set the Time Measure object. Off by default. 
       *
-       * @param time_measure
+       * @param time_measure 
       */
      void setTimeMeasure(TIME_MEASURE time_measure);

      /**
       * @brief Setting the time step for each iteration step. Time step must be
-       *        greater than zero.
+       *        greater than zero. Setting the timestep is required. 
       *
       * @param timestep Valid timestep greater than zero. 
       */
@ -127,7 +131,7 @@ class Simulation {

      /**
       * @brief Set the desired iterations to be calculated. A value greater
-       *        than zero must be specified here.
+       *        than zero must be specified here. Setting iterations is required. 
       *
       * @param iterations Number of iterations to be simulated.
       */
@ -164,6 +168,7 @@ class Simulation {
       */
      void printConcentrationsConsole();

+    // TODO move create CSVfile to TugUtils
      /**
       * @brief Creates a CSV file with a name containing the current simulation
       *        parameters. If the data name already exists, an additional counter is
--- a/src/BTCSv2.cpp
+++ b/src/BTCSv2.cpp
@ -1,7 +1,9 @@
 /**
- * @file BTCS.cpp
+ * @file BTCSv2.cpp
 * @brief Implementation of heterogenous BTCS (backward time-centered space) solution
- *        of diffusion equation in 1D and 2D space.
+ *        of diffusion equation in 1D and 2D space. Internally the alternating-direction
+ *        implicit (ADI) method is used. Version 2, because Version 1 was an 
+ *        implementation for the homogeneous BTCS solution. 
 * 
 */

@ -15,7 +17,7 @@ using namespace Eigen;


 // calculates coefficient for left boundary in constant case
-static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, double &sx) {
+static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, double sx) {
    double centerCoeff;
    double rightCoeff;

@ -28,7 +30,7 @@ static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int


 // calculates coefficient for left boundary in closed case
-static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, double &sx) {
+static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, double sx) {
    double centerCoeff;
    double rightCoeff;

@ -40,7 +42,7 @@ static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int &r


 // calculates coefficient for right boundary in constant case
-static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
+static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int rowIndex, int n, double sx) {
    double leftCoeff;
    double centerCoeff;

@ -53,7 +55,7 @@ static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int


 // calculates coefficient for right boundary in closed case
-static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
+static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int rowIndex, int n, double sx) {
    double leftCoeff;
    double centerCoeff;

@ -65,7 +67,7 @@ static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &


 // creates coefficient matrix for next time step from alphas in x-direction
-static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, int &numCols, int &rowIndex, double &sx) {
+static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, int numCols, int rowIndex, double sx) {

    // square matrix of column^2 dimension for the coefficients
    SparseMatrix<double> cm(numCols, numCols);
@ -119,7 +121,7 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryEl

 // calculates explicity concentration at top boundary in constant case
 static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrations, 
-                            MatrixXd &alpha, vector<BoundaryElement> &bcTop, int &rowIndex, int &i, double &sy) {
+                            MatrixXd &alpha, vector<BoundaryElement> &bcTop, int rowIndex, int i, double sy) {
    double c;

    c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
@ -138,7 +140,7 @@ static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrat

 // calculates explicit concentration at top boundary in closed case
 static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentrations, 
-                            MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
+                            MatrixXd &alpha, int rowIndex, int i, double sy) {
    double c;

    c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
@ -155,7 +157,7 @@ static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentratio

 // calculates explicit concentration at bottom boundary in constant case
 static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concentrations, 
-                            MatrixXd &alpha, vector<BoundaryElement> &bcBottom, int &rowIndex, int &i, double &sy) {
+                            MatrixXd &alpha, vector<BoundaryElement> &bcBottom, int rowIndex, int i, double sy) {
    double c;

    c = sy * alpha(rowIndex,i) * bcBottom[i].getValue()
@ -174,7 +176,7 @@ static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concent

 // calculates explicit concentration at bottom boundary in closed case
 static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentrations, 
-                            MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
+                            MatrixXd &alpha, int rowIndex, int i, double sy) {
    double c;

    c = (
@ -193,7 +195,7 @@ static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentra
 static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, 
                                        vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, 
                                        vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom, 
-                                        int &length, int &rowIndex, double &sx, double &sy) {
+                                        int length, int rowIndex, double sx, double sy) {

    VectorXd sv(length);
    int numRows = concentrations.rows();
@ -217,7 +219,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
    }

    // first row
-    if (rowIndex == 0) {
+    else if (rowIndex == 0) {
        for (int i = 0; i < length; i++) {
            type = bcTop[i].getType();
            if (type == BC_TYPE_CONSTANT) {
@ -231,7 +233,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
    }

    // last row
-    if (rowIndex == numRows-1) {
+    else if (rowIndex == numRows-1) {
        for (int i = 0; i < length; i++) {
            type = bcBottom[i].getType();
            if (type == BC_TYPE_CONSTANT) {
@ -316,7 +318,7 @@ static VectorXd ThomasAlgorithm(SparseMatrix<double> &A, VectorXd &b) {


 // BTCS solution for 1D grid 
-static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b)) {
+static void BTCS_1D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b)) {
    int length = grid.getLength();
    double sx = timestep / (grid.getDelta() * grid.getDelta());

@ -353,7 +355,7 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solve


 // BTCS solution for 2D grid
-static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b), int &numThreads) {
+static void BTCS_2D(Grid &grid, Boundary &bc, double timestep, VectorXd (*solverFunc) (SparseMatrix<double> &A, VectorXd &b), int numThreads) {
    int rowMax = grid.getRow();
    int colMax = grid.getCol();
    double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
@ -373,6 +375,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solve
    vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);

    MatrixXd concentrations = grid.getConcentrations();
+
    #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
    for (int i = 0; i < rowMax; i++) {
      
@ -394,7 +397,6 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solve
    alphaY.transposeInPlace();
    
    #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
-    
    for (int i = 0; i < colMax; i++) {
        // swap alphas, boundary conditions and sx/sy for column-wise calculation
        A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy);
@ -413,7 +415,7 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep, VectorXd (*solve


 // entry point for EigenLU solver; differentiate between 1D and 2D grid
-static void BTCS_LU(Grid &grid, Boundary &bc, double &timestep, int &numThreads) {
+static void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) {
    if (grid.getDim() == 1) {
        BTCS_1D(grid, bc, timestep, EigenLUAlgorithm);
    } else if (grid.getDim() == 2) {
@ -424,7 +426,7 @@ static void BTCS_LU(Grid &grid, Boundary &bc, double &timestep, int &numThreads)
 }

 // entry point for Thomas algorithm solver; differentiate 1D and 2D grid
-static void BTCS_Thomas(Grid &grid, Boundary &bc, double &timestep, int &numThreads) {
+static void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) {
    if (grid.getDim() == 1) {
        BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
    } else if (grid.getDim() == 2) {
--- a/src/Simulation.cpp
+++ b/src/Simulation.cpp
@ -232,6 +232,7 @@ void Simulation::run() {
    auto begin = std::chrono::high_resolution_clock::now();

    if (approach == FTCS_APPROACH) { // FTCS case
+
        for (int i = 0; i < iterations * innerIterations; i++) {
            if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
                printConcentrationsConsole();
@ -278,6 +279,10 @@ void Simulation::run() {
            }
        }

+    } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case
+
+        // TODO 
+
    }

    auto end = std::chrono::high_resolution_clock::now();
--- a/utils/output-visualization.html
+++ b/utils/output-visualization.html
@ -0,0 +1,543 @@
+<!DOCTYPE html>
+<html lang="en-US">
+<head>
+    <!-- <link rel="stylesheet" href="style.css"> -->
+    <style>
+        body {
+            background-color: white;
+        }
+
+        :root {
+            --max-width: 1440px;
+        }
+
+        .main-wrapper {
+            display: grid;
+            grid-template-areas: "sidebar main";
+            grid-template-columns: minmax(0,1fr) minmax(0,2.5fr);
+            grid-column-gap: 50px;
+
+            margin: 0 auto;
+            max-width: var(--max-width);
+        }
+
+        #state-info {
+            background-color: lightgray;
+            width: 120px;
+            text-align: center;
+            border: 2px solid black;
+            border-radius: 5px;
+            font-family: Arial, Helvetica, sans-serif;
+            font-size: 25px;
+            margin-bottom: 20px;
+        }
+    </style>
+</head>
+<body>
+<div id="root">
+<div class="main-wrapper">
+
+    <div class="sidebar-container">
+        <aside id="sidebar">
+        <div class="menu"> 
+            <form>
+                <label for="c_file">File (XTREME-format): </label>
+                <input type="file" id="c_file" name="c_file" accept="text/csv" />
+            </form>
+        </div>
+        <div class="info">
+            <p>
+                Use the keys <i>q</i> and <i>e</i> for increasing and decreasing the state by <b>1</b>.
+                The keys <i>a</i> and <i>d</i> for increasing and decresing the state by <b>10</b>.
+            </p>
+            <div id="state-info">
+                State: -
+            </div>
+            <div id="legend-info">
+
+            </div>
+        </div>
+        </aside>
+    </div>
+
+    <div class="content-container">
+        <main id="content" class="main-content">
+            
+        </main>
+    </div>
+
+</div>
+</div>
+</body>
+</html>
+
+<script type="module">
+
+    import * as d3 from "https://cdn.jsdelivr.net/npm/d3@7/+esm";
+
+    // layout information; constants
+    const width = document.getElementById('content').offsetWidth;
+    const height = document.getElementById('content').offsetWidth;
+    const svgMargin = { top: 10, right: 10, bottom: 10, left: 10 };
+    const boundary_gap = 1;
+    const gridWidth = width - svgMargin.left - svgMargin.right - 2*boundary_gap;
+    const gridHeight = height - svgMargin.top - svgMargin.bottom - 2*boundary_gap;
+
+    // simulation information
+    var state = 0; // iteration state
+    var iteration_count = 0;
+    var number_rows = 0;
+    var number_cols = 0;
+    var max_concentration = 0;
+    var data = [];
+    var leftb, rightb, topb, bottomb;
+
+    // test data; to be deleted later on
+    var test_data = [{"row": 0, "col": 0, "value": 1},
+                     {"row": 0, "col": 1, "value": 1},
+                     {"row": 0, "col": 2, "value": 1},
+                     {"row": 0, "col": 3, "value": 1},
+                     {"row": 1, "col": 0, "value": 1},
+                     {"row": 1, "col": 1, "value": 1},
+                     {"row": 1, "col": 2, "value": 1},
+                     {"row": 1, "col": 3, "value": 1},
+                     {"row": 2, "col": 0, "value": 1},
+                     {"row": 2, "col": 1, "value": 1},
+                     {"row": 2, "col": 2, "value": 1},
+                     {"row": 2, "col": 3, "value": 1},
+                     {"row": 3, "col": 0, "value": 1},
+                     {"row": 3, "col": 1, "value": 1},
+                     {"row": 3, "col": 2, "value": 1},
+                     {"row": 3, "col": 3, "value": 1}];
+
+
+    //////////////////////////////////////////
+    ////////
+    //////////////////////////////////////////
+    // Copyright 2021, Observable Inc.
+    // Released under the ISC license.
+    // https://observablehq.com/@d3/color-legend
+    function Legend(color, {
+    title,
+    tickSize = 6,
+    width = 320, 
+    height = 44 + tickSize,
+    marginTop = 18,
+    marginRight = 0,
+    marginBottom = 16 + tickSize,
+    marginLeft = 0,
+    ticks = width / 64,
+    tickFormat,
+    tickValues
+    } = {}) {
+
+    function ramp(color, n = 256) {
+        const canvas = document.createElement("canvas");
+        canvas.width = n;
+        canvas.height = 1;
+        const context = canvas.getContext("2d");
+        for (let i = 0; i < n; ++i) {
+        context.fillStyle = color(i / (n - 1));
+        context.fillRect(i, 0, 1, 1);
+        }
+        return canvas;
+    }
+
+    const svg = d3.create("svg")
+        .attr("width", width)
+        .attr("height", height)
+        .attr("viewBox", [0, 0, width, height])
+        .style("overflow", "visible")
+        .style("display", "block");
+
+    let tickAdjust = g => g.selectAll(".tick line").attr("y1", marginTop + marginBottom - height);
+    let x;
+
+    // Continuous
+    if (color.interpolate) {
+        const n = Math.min(color.domain().length, color.range().length);
+
+        x = color.copy().rangeRound(d3.quantize(d3.interpolate(marginLeft, width - marginRight), n));
+
+        svg.append("image")
+            .attr("x", marginLeft)
+            .attr("y", marginTop)
+            .attr("width", width - marginLeft - marginRight)
+            .attr("height", height - marginTop - marginBottom)
+            .attr("preserveAspectRatio", "none")
+            .attr("xlink:href", ramp(color.copy().domain(d3.quantize(d3.interpolate(0, 1), n))).toDataURL());
+    }
+
+    // Sequential
+    else if (color.interpolator) {
+        x = Object.assign(color.copy()
+            .interpolator(d3.interpolateRound(marginLeft, width - marginRight)),
+            {range() { return [marginLeft, width - marginRight]; }});
+
+        svg.append("image")
+            .attr("x", marginLeft)
+            .attr("y", marginTop)
+            .attr("width", width - marginLeft - marginRight)
+            .attr("height", height - marginTop - marginBottom)
+            .attr("preserveAspectRatio", "none")
+            .attr("xlink:href", ramp(color.interpolator()).toDataURL());
+
+        // scaleSequentialQuantile doesn’t implement ticks or tickFormat.
+        if (!x.ticks) {
+        if (tickValues === undefined) {
+            const n = Math.round(ticks + 1);
+            tickValues = d3.range(n).map(i => d3.quantile(color.domain(), i / (n - 1)));
+        }
+        if (typeof tickFormat !== "function") {
+            tickFormat = d3.format(tickFormat === undefined ? ",f" : tickFormat);
+        }
+        }
+    }
+
+    // Threshold
+    else if (color.invertExtent) {
+        const thresholds
+            = color.thresholds ? color.thresholds() // scaleQuantize
+            : color.quantiles ? color.quantiles() // scaleQuantile
+            : color.domain(); // scaleThreshold
+
+        const thresholdFormat
+            = tickFormat === undefined ? d => d
+            : typeof tickFormat === "string" ? d3.format(tickFormat)
+            : tickFormat;
+
+        x = d3.scaleLinear()
+            .domain([-1, color.range().length - 1])
+            .rangeRound([marginLeft, width - marginRight]);
+
+        svg.append("g")
+        .selectAll("rect")
+        .data(color.range())
+        .join("rect")
+            .attr("x", (d, i) => x(i - 1))
+            .attr("y", marginTop)
+            .attr("width", (d, i) => x(i) - x(i - 1))
+            .attr("height", height - marginTop - marginBottom)
+            .attr("fill", d => d);
+
+        tickValues = d3.range(thresholds.length);
+        tickFormat = i => thresholdFormat(thresholds[i], i);
+    }
+
+    // Ordinal
+    else {
+        x = d3.scaleBand()
+            .domain(color.domain())
+            .rangeRound([marginLeft, width - marginRight]);
+
+        svg.append("g")
+        .selectAll("rect")
+        .data(color.domain())
+        .join("rect")
+            .attr("x", x)
+            .attr("y", marginTop)
+            .attr("width", Math.max(0, x.bandwidth() - 1))
+            .attr("height", height - marginTop - marginBottom)
+            .attr("fill", color);
+
+        tickAdjust = () => {};
+    }
+
+    svg.append("g")
+        .attr("transform", `translate(0,${height - marginBottom})`)
+        .call(d3.axisBottom(x)
+            .ticks(ticks, typeof tickFormat === "string" ? tickFormat : undefined)
+            .tickFormat(typeof tickFormat === "function" ? tickFormat : undefined)
+            .tickSize(tickSize)
+            .tickValues(tickValues))
+        .call(tickAdjust)
+        .call(g => g.select(".domain").remove())
+        .call(g => g.append("text")
+            .attr("x", marginLeft)
+            .attr("y", marginTop + marginBottom - height - 6)
+            .attr("fill", "currentColor")
+            .attr("text-anchor", "start")
+            .attr("font-weight", "bold")
+            .attr("class", "title")
+            .text(title));
+
+    return svg.node();
+    }
+    //////////////////////////////////////////
+    ////////
+    //////////////////////////////////////////
+
+    
+
+    // max concentration
+    function calc_max_concentration(state) {
+        // var maxRow = data[state].map(function(row){ return Math.max.apply(Math, row); });
+        // var maxC = Math.max.apply(null, maxRow);
+        var maxC = Math.max(...data[state].map(d => d.value))
+        return maxC;
+    }
+
+    // scales
+    var col, row, color;
+    function create_n_update_scales() {
+        col = d3.scaleBand()
+            .range([0, gridWidth])
+            .domain([...Array(number_cols).keys()])
+            .padding(0.05);
+        row = d3.scaleBand()
+            .range([0, gridHeight])
+            .domain([...Array(number_rows).keys()])
+            .padding(0.05);
+        color = d3.scaleLinear()
+            .range(["#d0d0ec", "#7f0000"])
+            .domain([0, max_concentration]);
+    }
+
+    // grid
+    function draw_n_update_grid(grid_data) {
+        document.getElementById('state-info').innerHTML = `State: ${state}`;
+        // Grid rects
+        grid.selectAll()
+            .data(grid_data, (d,i) => d.row+':'+d.col)
+            .enter()
+            .append("rect")
+                .attr("x", (d,i) => col(d.col))
+                .attr("y", (d,i) => row(d.row))
+                .attr("width", col.bandwidth())
+                .attr("height", row.bandwidth())
+                .attr("fill", (d,i) => color(d.value))
+            .on("mouseover", mouseover)
+            .on("mousemove", mousemove)
+            .on("mouseleave", mouseleave);
+    }
+
+    // boundaries
+    function draw_boundaries() {
+        left_side.selectAll("rect")
+            .data(leftb)
+            .join("rect")
+            .attr("x", (d,i) => 0)
+            .attr("y", (d,i) => row(i))
+            .attr("width", svgMargin.left)
+            .attr("height", row.bandwidth())
+            .attr("fill", function(d,i) {
+                if (d == -1) {
+                    return "#000000";
+                } else {
+                    return color(d);
+                }
+            });
+        right_side.selectAll("rect")
+            .data(rightb)
+            .join("rect")
+            .attr("x", (d,i) => 0)
+            .attr("y", (d,i) => row(i))
+            .attr("width", svgMargin.right)
+            .attr("height", row.bandwidth())
+            .attr("fill", function(d,i) {
+                if (d == -1) {
+                    return "#000000";
+                } else {
+                    return color(d);
+                }
+            });
+        top_side.selectAll("rect")
+            .data(topb)
+            .join("rect")
+            .attr("x", (d,i) => col(i))
+            .attr("y", (d,i) => 0)
+            .attr("width", col.bandwidth())
+            .attr("height", svgMargin.top)
+            .attr("fill", function(d,i) {
+                if (d == -1) {
+                    return "#000000";
+                } else {
+                    return color(d);
+                }
+            });
+        bottom_side.selectAll("rect")
+            .data(bottomb)
+            .join("rect")
+            .attr("x", (d,i) => col(i))
+            .attr("y", (d,i) => 0)
+            .attr("width", col.bandwidth())
+            .attr("height", svgMargin.bottom)
+            .attr("fill", function(d,i) {
+                if (d == -1) {
+                    return "#000000";
+                } else {
+                    return color(d);
+                }
+            });
+    }
+
+    function draw_legend() {
+        d3.select("#legend-info").selectAll("svg").remove();
+
+        var legend = Legend(d3.scaleLinear([0, max_concentration], ["#d0d0ec", "#7f0000"]), {
+            title: "Concentration",
+            tickFormat: ".2f",
+            tickValues: [0, max_concentration/4, max_concentration/2, max_concentration/4*3, max_concentration]
+        });
+
+        d3.select("#legend-info")
+            .node()
+            .appendChild(legend);
+    }
+
+    // tooltip events
+    var mouseover = function(event, d) {
+        tooltip.style("opacity", 1);
+    }
+    var mousemove = function(event, d) {
+        const tooltip_text = d3.select('.tooltip-text');
+        const tooltip_box = d3.select('.tooltip-box');
+
+        tooltip_text.text(`Value: ${d.value}`);
+        var box_width = tooltip_text.node().getComputedTextLength();
+        tooltip_box.attr("width", box_width+10);
+
+        const [x, y] = d3.pointer(event);
+
+        tooltip
+          .attr('transform', `translate(${x+svgMargin.left+50}, ${y+svgMargin.top})`);
+    }
+    var mouseleave = function(event, d) {
+        tooltip.style("opacity", 0)
+    }
+
+
+    // key events for changing grid iteration state
+    addEventListener("keydown", (event) => {
+    if (event.isComposing || event.keyCode === 81) {
+        if (state > 0) {
+            state -= 1;
+            max_concentration = calc_max_concentration(state);
+            create_n_update_scales();
+            draw_legend();
+            draw_n_update_grid(data[state]);
+        }
+    }
+    if (event.isComposing || event.keyCode === 69) {
+        if (state < iteration_count-1) {
+            state += 1;
+            max_concentration = calc_max_concentration(state);
+            create_n_update_scales();
+            draw_legend();
+            draw_n_update_grid(data[state]);
+        }
+    }
+    if (event.isComposing || event.keyCode === 65) {
+        if (state > 9) {
+            state -= 10;
+        } else {
+            state = 0;
+        }
+        max_concentration = calc_max_concentration(state);
+        create_n_update_scales();
+        draw_legend();
+        draw_n_update_grid(data[state]);
+    }
+    if (event.isComposing || event.keyCode === 68) {
+        if (state < iteration_count-10) {
+            state += 10;
+        } else {
+            state = iteration_count-1;
+        }
+        max_concentration = calc_max_concentration(state);
+        create_n_update_scales();
+        draw_legend();
+        draw_n_update_grid(data[state]);
+    }
+    });
+
+
+    // read data from file when selecting new file
+    async function create_data_struct(csv_text) {
+        var row_data = d3.dsvFormat(" ").parseRows(csv_text);
+
+        leftb = row_data[0];
+        rightb = row_data[1];
+        topb = row_data[2];
+        bottomb = row_data[3];
+
+        number_rows = row_data[0].length;
+        number_cols = row_data[2].length;
+        iteration_count = (row_data.length - 6) / (number_rows + 2);
+        console.log(iteration_count);
+        var iteration = []; // loop temp variable
+        var entry = {}; // loop temp variable
+        for (let i = 0; i < iteration_count; i++) {
+            iteration = []; // reset temp variable
+            for (let j = 0; j < number_rows; j++) {
+                for (let k = 0; k < number_cols; k++) {
+                    entry = {}; // reset temp variable
+                    entry.row = j;
+                    entry.col = k;
+                    entry.value = row_data[6 + i*(number_rows+2) + j][k];
+                    iteration.push(entry);
+                }
+            }
+            data.push(iteration);
+        }
+    }
+    function reader(file, callback) {
+        const fr = new FileReader();
+        fr.onload = () => callback(null, fr.result);
+        fr.onerror = (err) => callback(err);
+        fr.readAsText(file)
+    }
+    document.getElementById("c_file").addEventListener("change", (evt) => {
+        if (!evt.target.files) { // no file, do nothing
+            return;
+        }
+        reader(evt.target.files[0], (err, res) => {
+            create_data_struct(res);
+            state = 0;
+            max_concentration = calc_max_concentration(state);
+            create_n_update_scales();
+            draw_legend();
+            draw_boundaries();
+            draw_n_update_grid(data[0]);
+        });
+    });
+
+
+    const svg = d3.create("svg")
+        .attr("width", width)
+        .attr("height", height);
+
+    var left_side = svg.append("g")
+        .attr("class", "boundary")
+        .attr("transform", `translate(0, ${svgMargin.top+boundary_gap})`);
+    var right_side = svg.append("g")
+        .attr("class", "boundary")
+        .attr("transform", `translate(${width-svgMargin.right}, ${svgMargin.top+boundary_gap})`);
+    var top_side = svg.append("g")
+        .attr("class", "boundary")
+        .attr("transform", `translate(${svgMargin.left+boundary_gap}, 0)`);
+    var bottom_side = svg.append("g")
+        .attr("class", "boundary")
+        .attr("transform", `translate(${svgMargin.left+boundary_gap}, ${height-svgMargin.bottom})`);
+
+    var grid = svg.append("g")
+        .attr("class", "grid")
+        .attr("transform", `translate(${svgMargin.left + boundary_gap}, ${svgMargin.top + boundary_gap})`);
+        
+    // tooltip 
+    var tooltip = svg.append("g")
+        .attr("class", "tooltip-area")
+        .style("opacity", 0);
+    var tooltip_box = tooltip.append("rect")
+        .attr("class", "tooltip-box")
+        .attr("height", 25)
+        .attr("transform", `translate(-5,-18)`)
+        .attr("fill", "white")
+        .attr("stroke", "black");
+    tooltip.append("text")
+        .attr("class", "tooltip-text")
+        .text("Value: ");
+    
+    content.append(svg.node());
+
+</script>