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

.gitlab-ci.yml

@@ -1,4 +1,4 @@
-image: sobc/gitlab-ci
+image: git.gfz-potsdam.de:5000/naaice/tug:ci
 
 stages:
     - build
@@ -6,8 +6,6 @@ stages:
     - static_analyze
 
 build_release:
-  before_script:
-    - apt-get update && apt-get install -y libeigen3-dev git
   stage: build
   artifacts:
       paths:
@@ -24,11 +22,11 @@ test:
         - ./build/test/testTug
 
 lint:
-  before_script:
-    - apt-get update && apt-get install -y libeigen3-dev
   stage: static_analyze
+  before_script:
+    - apk add clang-extra-tools
   script:
     - mkdir lint && cd lint
-    - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*,readability-*, modernize-*" -DTUG_ENABLE_TESTING=OFF ..
+    - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*, modernize-*" -DTUG_ENABLE_TESTING=OFF ..
     - make tug
   when: manual

doc/ADI_scheme.org

@@ -1,7 +1,7 @@
-#+TITLE: Numerical solution of diffusion equation in 2D with ADI Scheme
+#+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{amsmath, systeme}
+#+LATEX_HEADER: \usepackage{amsmath, systeme, cancel, xcolor}
 #+OPTIONS: toc:nil
 
 
@@ -308,8 +308,8 @@ form:
 
 * Heterogeneous diffusion
 
-If the diffusion coefficient $\alpha$ is spatially variable, [[eqn:1]] can
+If the diffusion coefficient $\alpha$ is spatially variable, equation
-be rewritten:
+[[eqn:1]] can be rewritten as:
 
 #+NAME: eqn:hetdiff
 \begin{align}
@@ -391,8 +391,8 @@ concentrations.
 ** Direct discretization
 
 As noted in literature (LeVeque and Numerical Recipes) a better way is
-to discretize directly the physical problem ([[eqn:hetdiff]]) at points
+to discretize directly the physical problem (eq. [[eqn:hetdiff]]) at
-halfway between grid points:
+points halfway between grid points:
 
 \begin{align*}
 \begin{cases}
@@ -401,15 +401,18 @@ halfway between grid points:
 \end{cases}
 \end{align*}
 
-\noindent A further differentiation gives us the centered
+\noindent A further differentiation gives us the spatially centered
 approximation of $\frac{\partial}{\partial x} \left(\alpha(x)
 \frac{\partial C }{\partial x}\right)$:
 
-\begin{align*}
+#+NAME: eqn:CS_het
+\begin{equation}
+\begin{aligned}
 \frac{\partial}{\partial x} \left(\alpha(x)
 \frac{\partial C }{\partial x}\right)(x_i) & \simeq \frac{1}{\Delta x}\left[\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) - \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \right]\\
 &\displaystyle =\frac{1}{\Delta x^2} \left[ \alpha_{i+1/2}C_{i+1} - (\alpha_{i+1/2}+\alpha_{i-1/2}) C_{i} + \alpha_{i-1/2}C_{i-1}\right]
-\end{align*}
+\end{aligned}
+\end{equation}
 
 \noindent The ADI scheme with this approach becomes:
 
@@ -439,3 +442,98 @@ explicit terms, the two sweeps become:
 \end{aligned}
 \right.
 \end{equation}
+
+\bigskip
+
+\noindent The "interblock" diffusion coefficients $\alpha_{i+1/2,j}$
+can be arithmetic mean:
+
+\[
+\displaystyle \alpha_{i+1/2, j} = \displaystyle \frac{\alpha_{i+1, j} + \alpha_{i, j}}{2}
+\]
+
+\noindent or the harmonic mean:
+
+\[
+\displaystyle \alpha_{i+1/2, j} = \displaystyle \frac{2}{\frac{1}{\alpha_{i+1, j}} + \frac{1}{\alpha_{i, j}}}
+\]
+
+
+
+\pagebreak
+
+* Explicit scheme for 2D heterogeneous diffusion
+
+A classical explicit FTCS scheme (forward in time, central in space)
+for 2D heterogeneous diffusion can be expressed simply leveraging the
+discretization of equation [[eqn:CS_het]]:
+
+#+NAME: eqn:2DHeterFTCS
+\begin{equation}
+\begin{aligned}
+\frac{C_{i,j}^{t+1} - C_{i,j}^{t}}{\Delta t} = & \frac{1}{\Delta x^2} \left[ \alpha^x_{i+1/2, j}C^t_{i+1, j} - (\alpha^x_{i+1/2, j} + \alpha^x_{i-1/2, j}) C^t_{i,j} + \alpha^x_{i-1/2,j}C^t_{i-1,j}\right] + \\
+                                               & \frac{1}{\Delta y^2} \left[ \alpha^y_{i, j+1/2}C^t_{i, j+1} - (\alpha^y_{i, j+1/2} + \alpha^y_{i, j-1/2}) C^t_{i,j} + \alpha^y_{i,j-1/2}C^t_{i,j-1}\right]
+\end{aligned}
+\end{equation}
+\noindent where in the RHS only the known concentrations at time $t$
+appear. Rearranging the terms, we get:
+
+#+NAME: eqn:2DHeterFTCS_final
+\begin{equation}
+\begin{aligned}
+C_{i,j}^{t+1} = & C^t_{i,j} +\\
+                & \frac{\Delta t}{\Delta x^2} \left[ \alpha^x_{i+1/2, j}C^t_{i+1, j} - (\alpha^x_{i+1/2, j} + \alpha^x_{i-1/2, j}) C^t_{i,j} + \alpha^x_{i-1/2,j}C^t_{i-1,j}\right] + \\
+                & \frac{\Delta t}{\Delta y^2} \left[ \alpha^y_{i, j+1/2}C^t_{i, j+1} - (\alpha^y_{i, j+1/2} + \alpha^y_{i, j-1/2}) C^t_{i,j} + \alpha^y_{i,j-1/2}C^t_{i,j-1}\right]
+\end{aligned}
+\end{equation}
+
+The Courant-Friedrichs-Lewy stability criterion for this scheme reads:
+#+NAME: eqn:CFL2DFTCS
+\begin{equation}
+\Delta t \leq \frac{(\Delta x^2, \Delta y^2)}{2 \max(\alpha_{i,j})}
+\end{equation}
+
+** Boundary conditions
+
+In analogy to the treatment of the 1D homogeneous FTCS scheme (cfr
+section 1), we need to differentiate the domain boundaries ($i=0$ and
+$i=n_x$; the same applies to $j$ of course) accounting for the
+discrepancy in the discretization.
+
+For the zero-th (left) cell, whose center is at $x=dx/2$, we can
+evaluate the left gradient with the left boundary using such distance,
+calling $l$ the numerical value of a constant boundary condition,
+equation [[eqn:CS_het]] becomes:
+
+#+NAME: eqn:2D_FTCS_left
+\begin{equation}
+\begin{aligned}
+\frac{\partial}{\partial x} \left(\alpha(x)
+\frac{\partial C }{\partial x}\right)(x_0) & \simeq \frac{1}{\Delta x}\left[\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) - \alpha_{i} \left( \frac{C_{i} - l }{\frac{\Delta x}{2}} \right) \right]\\
+&\displaystyle =\frac{1}{\Delta x^2} \left[ \alpha_{i+1/2}C_{i+1} - (\alpha_{i+1/2}+ 2\alpha_i) C_{i} + 2 \alpha_{i}\cdot l\right]
+\end{aligned}
+\end{equation}
+
+\noindent Similarly, for $i=n_x$,
+
+#+NAME: eqn:2D_FTCS_right
+\begin{equation}
+\begin{aligned}
+\frac{\partial}{\partial x} \left(\alpha(x)
+\frac{\partial C }{\partial x}\right)(x_n) & \simeq \frac{1}{\Delta x}\left[\alpha_{i} \left( \frac{r -C_{i}}{\frac{\Delta x}{2}} \right) - \alpha_{i-1/2} \left( \frac{C_{i} - C_{i-1}} {\Delta x} \right) \right]\\
+&\displaystyle =\frac{1}{\Delta x^2} \left[ 2 \alpha_{i} r - (\alpha_{i+1/2}+ 2\alpha_i) C_{i} + \alpha_{i-1/2}\cdot C_{i-1}\right]
+\end{aligned}
+\end{equation}
+
+If on the right boundary we have *closed* or Neumann condition, the
+left derivative becomes zero and we are left with:
+
+#+NAME: eqn:2D_FTCS_rightclosed
+\begin{equation}
+\begin{aligned}
+\frac{\partial}{\partial x} \left(\alpha(x)
+\frac{\partial C }{\partial x}\right)(x_n) & \simeq \frac{1}{\Delta x}\left[\cancel{\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right)} - \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \right]\\
+&\displaystyle =\frac{\alpha_{i-1/2}}{\Delta x^2} (C_{i-1} - C_i)
+\end{aligned}
+\end{equation}
+

examples/CMakeLists.txt

@@ -2,9 +2,11 @@ add_executable(first_example first_example.cpp)
 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)
 
 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(FTCS_2D_proto_example easy_profiler) 

examples/FTCS_1D_proto_example.cpp

@@ -0,0 +1,44 @@
+#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,20);
+    grid.setConcentrations(concentrations);
+
+
+    // ******************
+    // **** BOUNDARY ****
+    // ******************
+
+    // create a boundary with constant values
+    Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
+
+
+    // ************************
+    // **** SIMULATION ENV ****
+    // ************************
+
+    // set up a simulation environment
+    Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // grid,boundary,simulation-approach
+
+    // (optional) set the timestep of the simulation
+    simulation.setTimestep(0.1); // timestep
+
+    // (optional) set the number of iterations
+    simulation.setIterations(100);
+
+    // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
+    simulation.setOutputCSV(CSV_OUTPUT_OFF);
+    
+    // **** RUN SIMULATION ****
+
+    // run the simulation
+    simulation.run();
+}

examples/FTCS_2D_proto_example.cpp

@@ -49,12 +49,12 @@ int main(int argc, char *argv[]) {
     // (optional) set boundary condition values for one side, e.g.:
     // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value
     // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values
-    VectorXd bc_zero_values = VectorXd::Constant(20,0);
+    // VectorXd bc_zero_values = VectorXd::Constant(20,0);
-    bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values);
+    // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values);
-    bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values);
+    // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values);
-    VectorXd bc_front_values = VectorXd::Constant(20,2000);
+    // VectorXd bc_front_values = VectorXd::Constant(20,2000);
-    bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values);
+    // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values);
-    bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values);
+    // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values);
 
 
     // ************************

include/tug/Boundary.hpp

@@ -4,6 +4,7 @@
 #include <Eigen/Core>
 #include "Grid.hpp"
 
+using namespace std;
 using namespace Eigen;
 
 enum BC_TYPE {
@@ -18,6 +19,45 @@ enum BC_SIDE {
     BC_SIDE_BOTTOM
 };
 
+/*******************************************
+class WallElement {
+    public:
+        WallElement() {
+            this->type = BC_TYPE_CLOSED;
+            this->value = 0;
+        }
+
+        WallElement(double value) {
+            this->type = BC_TYPE_CONSTANT;
+            this->value = value;
+        }
+
+        BC_TYPE getType() {
+            return this->type;
+        }
+
+        double getValue() {
+            return this->value;
+        }
+
+    private:
+        BC_TYPE type;
+        double value;
+};
+
+class BoundaryWall {
+    public: 
+        BoundaryWall(int length) {
+            // create array with length many wall elements
+        }
+
+    private:
+        BC_SIDE side;
+        int length;
+        vector<WallElement> wall;
+};
+***********************/
+
 class Boundary {
     public:
 
@@ -55,6 +95,10 @@ class Boundary {
 
     private:
         Grid grid;
+
+        // need a way to save the bc type and value for each single 'boundary cell'
+        // perhaps an array for each side with structs containing the bc type as well as a value
+        // or another object that contains one boundary side 
         BC_TYPE type;
         VectorXd left, right, top, bottom;
 };

proto/FTCS.ipynb

@@ -458,7 +458,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.7"
+   "version": "3.11.4"
   },
   "orig_nbformat": 4
  },

scripts/Adi2D_Reference.R

@@ -1,4 +1,4 @@
-## Time-stamp: "Last modified 2023-01-05 17:52:55 delucia"
+## Time-stamp: "Last modified 2023-05-11 17:31:41 delucia"
 
 ## Brutal implementation of 2D ADI scheme
 ## Square NxN grid with dx=dy=1
@@ -272,8 +272,8 @@ ADIHetDir <- function(field, dt, iter, alpha) {
     
     for (it in seq(1, iter)) {
         for (i in seq(2, ny-1)) {
-            Aij <- cbind(harm(alpha[i,], alpha[i-1,]), harm(alpha[i,], alpha[i+1,]))
+            Aij <- cbind(colMeans(rbind(alpha[i,], alpha[i-1,])), colMeans(rbind(alpha[i,], alpha[i+1,])))
-            Bij <- cbind(harm(alpha[,i], alpha[,i-1]), harm(alpha[,i], alpha[,i+1]))
+            Bij <- cbind(rowMeans(cbind(alpha[,i], alpha[,i-1])), rowMeans(cbind(alpha[,i], alpha[,i+1])))
             tmpX[i,] <- SweepByRowHetDir(i, res, dt=dt, Aij, Bij)
         }
         resY <- t(tmpX)
@@ -321,22 +321,27 @@ SweepByRowHetDir <- function(i, field, dt, Aij, Bij) {
 ## adi2 <- ADI(n=51, dt=10, iter=200, alpha=1E-3)
 ## ref2 <- DoRef(n=51, alpha=1E-3, dt=10, iter=200)
 
-n <- 5
+n <- 51
 field <- matrix(0, n, n)
 alphas <- matrix(1E-5*runif(n*n, 1,2), n, n)
 
-alphas1 <- matrix(3E-5, n, 25)
-alphas2 <- matrix(1E-5, n, 26)
 
-alphas <- cbind(alphas1, alphas2)
+## dim(field)
+## dim(alphas)
+## all.equal(dim(field), dim(alphas))
+
+## alphas1 <- matrix(3E-5, n, 25)
+## alphas2 <- matrix(1E-5, n, 26)
+
+## alphas <- cbind(alphas1, alphas2)
 
 ## for (i in seq(1,nrow(alphas)))
 ##     alphas[i,] <- seq(1E-7,1E-3, length=n)
 
 #diag(alphas) <- rep(1E-2, n)
 
-adih  <- ADIHetDir(field=field, dt=10, iter=200, alpha=alphas)
+adih  <- ADIHetDir(field=field, dt=20, iter=500, alpha=alphas)
-adi2  <- ADI(n=n, dt=10, iter=200, alpha=1E-5)
+adi2  <- ADI(n=n, dt=20, iter=500, alpha=1E-5)
 
 
 par(mfrow=c(1,3))
@@ -347,6 +352,17 @@ plot(adih[[length(adih)]], adi2[[length(adi2)]], pch=4, log="xy")
 abline(0,1)
 
 
+cchet <- lapply(adih, round, digits=6)
+cchom <- lapply(adi2, round, digits=6)
+
+plot(cchet[[length(cchet)]], cchom[[length(cchom)]], pch=4, log="xy", xlim=c(1e-6,1), ylim=c(1e-6,1))
+abline(0,1)
+
+cchet[[500]]
+
+str(adih)
+
+
 sapply(adih, sum)
 sapply(adi2, sum)
 
@@ -360,3 +376,108 @@ image(adih[[length(adih)]])
 points(0.5,0.5, col="red",pch=4)
 
 options(width=110)
+
+
+FTCS_2D <- function(field, dt, iter, alpha) {
+
+    if (!all.equal(dim(field), dim(alpha)))
+        stop("field and alpha are not matrix")
+    
+    ## now both field and alpha must be nx*ny matrices
+    nx <- ncol(field)
+    ny <- nrow(field)
+    dx <- dy <- 1
+
+    ## find out the center of the grid to apply conc=1
+    cenx <- ceiling(nx/2)
+    ceny <- ceiling(ny/2)
+    field[cenx, ceny] <- 1
+
+    ## prepare containers for computations and outputs
+    tmp <- res <- field
+
+    cflt <- 1/max(alpha)/4
+    cat(":: CFL allowable time step: ", cflt,"\n")
+
+    ## inner iterations
+    inner <- floor(dt/cflt)
+    if (inner == 0) {
+        ## dt < cflt, no inner iterations
+        inner <- 1
+        tsteps <- dt
+        cat(":: No inner iter. required\n")
+    } else {
+        tsteps <- c(rep(cflt, inner), dt-inner*cflt)
+        cat(":: Number of inner iter. required: ", inner,"\n")
+    }
+    
+    
+    out <- vector(mode="list", length=iter)
+
+    for (it in seq(1, iter)) {
+        cat(":: outer it: ", it)
+        
+        for (innerit in seq_len(inner)) {
+            for (i in seq(2, ny-1)) {
+                for (j in seq(2, nx-1)) {
+                    ## tmp[i,j] <- res[i,j] +
+                    ## + tsteps[innerit]/dx/dx * (res[i+1,j]*mean(alpha[i+1,j],alpha[i,j]) -
+                    ##                            res[i,j]  *(mean(alpha[i+1,j],alpha[i,j])+mean(alpha[i-1,j],alpha[i,j])) +
+                    ##                            res[i-1,j]*mean(alpha[i-1,j],alpha[i,j])) +
+                    ## + tsteps[innerit]/dy/dy * (res[i,j+1]*mean(alpha[i,j+1],alpha[i,j]) -
+                    ##                            res[i,j]  *(mean(alpha[i,j+1],alpha[i,j])+mean(alpha[i,j-1],alpha[i,j])) +
+                    ##                            res[i,j-1]*mean(alpha[i,j-1],alpha[i,j]))
+                    tmp[i,j] <- res[i,j] +
+                        + tsteps[innerit]/dx/dx * ((res[i+1,j]-res[i,j]) * mean(alpha[i+1,j],alpha[i,j]) -
+                                                   (res[i,j]-res[i-1,j]) * mean(alpha[i-1,j],alpha[i,j])) +
+                        + tsteps[innerit]/dx/dx * ((res[i,j+1]-res[i,j]) * mean(alpha[i,j+1],alpha[i,j]) -
+                                                   (res[i,j]-res[i,j-1]) * mean(alpha[i,j-1],alpha[i,j])) 
+                }
+            }
+            ## swap back tmp to res for the next inner iteration 
+            res <- tmp
+        } 
+        cat("- done\n")
+        ## at end of inner it we store 
+        out[[it]] <- res
+    }
+    
+    return(out)
+}
+
+## testing that FTCS with homog alphas reverts to ADI/Reference sim
+n <- 51
+field <- matrix(0, n, n)
+alphas <- matrix(1E-3, n, n)
+
+adi2 <- ADI(n=51, dt=100, iter=20, alpha=1E-3)
+
+ref  <-  DoRef(n=51, alpha=1E-3, dt=100, iter=20)
+
+adihet  <- ADIHetDir(field=field, dt=100, iter=20, alpha=alphas)
+
+ftcsh <- FTCS_2D(field=field, dt=100, iter=20, alpha=alphas)
+
+
+par(mfrow=c(2,4))
+image(ref, main="Reference ODE.2D")
+points(0.5,0.5, col="red",pch=4)
+image(ftcsh[[length(ftcsh)]], main="FTCS 2D")
+points(0.5,0.5, col="red",pch=4)
+image(adihet[[length(adihet)]], main="ADI Heter.")
+points(0.5,0.5, col="red",pch=4)
+image(adi2[[length(adi2)]], main="ADI Homog.", col=terrain.colors(12))
+points(0.5,0.5, col="red",pch=4)
+plot(ftcsh[[length(ftcsh)]], ref, pch=4, log="xy", xlim=c(1E-16, 1), ylim=c(1E-16, 1),
+     main = "FTCS_2D vs ref", xlab="FTCS 2D", ylab="Reference")
+abline(0,1)
+plot(ftcsh[[length(ftcsh)]], adihet[[length(adihet)]], pch=4, log="xy", xlim=c(1E-16, 1), ylim=c(1E-16, 1),
+     main = "FTCS_2D vs ADI Het", xlab="FTCS 2D", ylab="ADI 2D Heter.")
+abline(0,1)
+plot(ftcsh[[length(ftcsh)]], adi2[[length(adi2)]], pch=4, log="xy", xlim=c(1E-16, 1), ylim=c(1E-16, 1),
+     main = "FTCS_2D vs ADI Hom", xlab="FTCS 2D", ylab="ADI 2D Hom.")
+abline(0,1)
+plot(adihet[[length(adihet)]], adi2[[length(adi2)]], pch=4, log="xy", xlim=c(1E-16, 1), ylim=c(1E-16, 1),
+     main = "ADI Het vs ADI Hom", xlab="ADI Het", ylab="ADI 2D Hom.")
+abline(0,1)
+

src/FTCS.cpp

@@ -4,7 +4,7 @@
 
 using namespace std;
 
-double calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = false) {
+double calcAlphaIntercell(double alpha1, double alpha2, bool useHarmonic = false) {
     if (useHarmonic) {
         return 2 / ((1/alpha1) + (1/alpha2));
     } else {
@@ -12,7 +12,168 @@ double calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = fal
     }
 }
 
-MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
+
+double calcHorizontalChange(Grid grid, int row, int col) {
+
+    double result = 
+        calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+            * grid.getConcentrations()(row,col+1)
+        - (
+            calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
+            + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+            )
+            * grid.getConcentrations()(row,col)
+        + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+            * grid.getConcentrations()(row,col-1);
+
+    return result;
+}
+
+
+double calcVerticalChange(Grid grid, int row, int col) {
+    
+    double result =    
+        calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+            * grid.getConcentrations()(row+1,col)
+        - (
+            calcAlphaIntercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+            + calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
+            )
+            * grid.getConcentrations()(row,col)
+        + calcAlphaIntercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
+            * grid.getConcentrations()(row-1,col);
+
+    return result;
+}
+
+
+double calcHorizontalChangeLeftBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
+
+    double result = 
+        calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+            * grid.getConcentrations()(row,col+1)
+        - (
+            calcAlphaIntercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+            + 2 * grid.getAlphaX()(row,col)
+            ) 
+            * grid.getConcentrations()(row,col) 
+        + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row);
+
+    return result;
+}
+
+
+double calcHorizontalChangeLeftBoundaryClosed() {
+    return 0;
+}
+
+
+double calcHorizontalChangeRightBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
+
+    double result = 
+        2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+        - (
+            calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+            + 2 * grid.getAlphaX()(row,col)
+            ) 
+            * grid.getConcentrations()(row,col)
+        + calcAlphaIntercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+            * grid.getConcentrations()(row,col-1);
+
+    return result;
+}
+
+
+double calcHorizontalChangeRightBoundaryClosed() {
+    return 0;
+}
+
+
+double calcVerticalChangeTopBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
+    
+    double result = 
+        calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) 
+            * grid.getConcentrations()(row+1,col)
+        - (
+            calcAlphaIntercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) 
+            + 2 * grid.getAlphaY()(row, col)
+            ) 
+            * grid.getConcentrations()(row, col)
+        + 2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col);
+
+    return result;
+}
+
+
+double calcVerticalChangeTopBoundaryClosed() {
+    return 0;
+}
+
+
+double calcVerticalChangeBottomBoundaryConstant(Grid grid, Boundary bc, int row, int col) {
+
+    double result = 
+        2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
+        - (
+            calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) 
+            + 2 * grid.getAlphaY()(row, col)
+            ) 
+            * grid.getConcentrations()(row, col)
+        + calcAlphaIntercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) 
+            * grid.getConcentrations()(row-1,col);
+
+    return result;
+}
+
+
+double calcVerticalChangeBottomBoundaryClosed() {
+    return 0;
+}
+
+
+MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
+    int colMax = grid.getCol();
+    double deltaCol = grid.getDeltaCol();
+
+    MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0);
+
+    // only one row in 1D case
+    int row = 0;
+
+    // inner cells
+    for (int col = 1; col < colMax-1; col++) {
+        concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+            + timestep / (deltaCol*deltaCol)
+                * (
+                    calcHorizontalChange(grid, row, col)
+                )
+            ;
+    }
+
+    // left boundary
+    int col = 0;
+    concentrations_t1(row, col) = grid.getConcentrations()(row,col)
+            + timestep / (deltaCol*deltaCol) 
+                * (
+                    calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
+                )
+            ;
+
+
+    // right boundary
+    col = colMax-1;
+    concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+            + timestep / (deltaCol*deltaCol) 
+                * (
+                    calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
+                )
+            ;
+
+    return concentrations_t1;
+}
+
+
+MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
     int rowMax = grid.getRow();
     int colMax = grid.getCol();
     double deltaRow = grid.getDeltaRow();
@@ -20,32 +181,20 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
 
     // Matrix with concentrations at time t+1
     // TODO profiler / only use 2 matrices
-    MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 1);
+    MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0);
 
     // inner cells
-    // (should have 7 calls to current concentration)
+    // these do not depend on the boundary condition type
     for (int row = 1; row < rowMax-1; row++) {
         for (int col = 1; col < colMax-1; col++) {
             concentrations_t1(row, col) = grid.getConcentrations()(row, col) 
                 + timestep / (deltaRow*deltaRow) 
                     * (
-                    calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+                        calcVerticalChange(grid, row, col)
-                    * grid.getConcentrations()(row+1,col)
-                    - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
-                    + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
-                    * grid.getConcentrations()(row,col)
-                    + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
-                    * grid.getConcentrations()(row-1,col)
                     )
                 + timestep / (deltaCol*deltaCol) 
                     * (
-                    calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+                        calcHorizontalChange(grid, row, col)
-                    * grid.getConcentrations()(row,col+1)
-                    - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
-                    + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)))
-                    * grid.getConcentrations()(row,col)
-                    + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                    * grid.getConcentrations()(row,col-1)
                     )
                 ;
         }
@@ -53,184 +202,119 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
 
     // boundary conditions
     // left without corners / looping over rows
-    // (should have 6 calls to current concentration)
     int col = 0;
     for (int row = 1; row < rowMax-1; row++) {
         concentrations_t1(row, col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
-                * (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+                * (
-                * grid.getConcentrations()(row,col+1)
+                    calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
-                - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+                )
-                + 2 * grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col) 
-                + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row))
             + timestep / (deltaRow*deltaRow)
-                * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+                * (
-                * grid.getConcentrations()(row+1,col)
+                    calcVerticalChange(grid, row, col)
-                - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+                )
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
+            ;
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
-                * grid.getConcentrations()(row-1,col));
     }
 
     // right without corners / looping over columns
-    // (should have 6 calls to current concentration)
     col = colMax-1;
     for (int row = 1; row < rowMax-1; row++) {
         concentrations_t1(row,col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
-                * (2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+                * (
-                - (calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+                    calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
-                + 2 * grid.getAlphaX()(row,col)) 
+                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                * grid.getConcentrations()(row,col-1))
             + timestep / (deltaRow*deltaRow)
-                * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+                * (
-                * grid.getConcentrations()(row+1,col)
+                    calcVerticalChange(grid, row, col)
-                - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+                )
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
+            ;
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
-                * grid.getConcentrations()(row-1,col));
     }
 
 
     // top without corners / looping over cols
-    // (should have 6 calls to current concentration)
     int row = 0;
     for (int col=1; col<colMax-1;col++){
         concentrations_t1(row, col) = grid.getConcentrations()(row, col)
-            + timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) * grid.getConcentrations()(1,col)
+            + timestep / (deltaRow*deltaRow) 
-                - (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) + 2 * grid.getAlphaY()(0, col)) * grid.getConcentrations()(0, col)
+                * (
-                + 2 * grid.getAlphaY()(0, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col))
+                    calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
-            + timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col+1)
+                )
-                - (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) + calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col))) * grid.getConcentrations()(0, col)
+            + timestep / (deltaCol*deltaCol) 
-                + calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col-1));
+                * (
+                    calcHorizontalChange(grid, row, col)
+                )
+            ;
     }
 
     // bottom without corners / looping over cols
-    // (should have 6 calls to current concentration)
     row = rowMax-1;
     for(int col=1; col<colMax-1;col++){
         concentrations_t1(row, col) = grid.getConcentrations()(row, col)
-            + timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
+            + timestep / (deltaRow*deltaRow) 
-                - (calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) + 2 * grid.getAlphaY()(row, col)) * grid.getConcentrations()(row, col)
+                * (
-                + calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) * grid.getConcentrations()(row-1,col))
+                    calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
-            + timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) * grid.getConcentrations()(row, col+1)
+                )
-                - (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) + calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col))) * grid.getConcentrations()(row, col)
+            + timestep / (deltaCol*deltaCol) 
-                + calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col)) * grid.getConcentrations()(row, col-1));
+                * (
+                    calcHorizontalChange(grid, row, col)
+                )
+            ;
     }
 
     // corner top left
-    // (should have 5 calls to current concentration)
     row = 0;
     col = 0;
     concentrations_t1(row,col) = grid.getConcentrations()(row,col)
         + timestep/(deltaCol*deltaCol)
             * (
-                calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col+1)
+                calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
-                - ( 
-                    calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
-                    + 2 * grid.getAlphaX()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) * grid.getConcentrations()(row+1,col)
+                calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
-                    + 2 * grid.getAlphaY()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + 2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col)
             )
         ;
 
     // corner top right
-    // (should have 5 calls to current concentration)
     row = 0;
     col = colMax-1;
     concentrations_t1(row,col) = grid.getConcentrations()(row,col) 
         + timestep/(deltaCol*deltaCol)
             * (
-                2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+                calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                    + 2 * grid.getAlphaX()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                * grid.getConcentrations()(row,col-1)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) * grid.getConcentrations()(row+1,col)
+                calcVerticalChangeTopBoundaryConstant(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
-                    + 2 * grid.getAlphaY()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + 2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col)
             )
         ;
 
     // corner bottom left
-    // (should have 5 calls to current concentration)
     row = rowMax-1;
     col = 0;
     concentrations_t1(row,col) = grid.getConcentrations()(row,col)
         + timestep/(deltaCol*deltaCol)
             * (
-                calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col+1)
+                calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col)
-                - ( 
-                    calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
-                    + 2 * grid.getAlphaX()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
+                calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaY()(row,col), grid.getAlphaY()(row-1,col))
-                    + 2 * grid.getAlphaY()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaY()(row,col), grid.getAlphaY()(row-1,col))
-                * grid.getConcentrations()(row-1,col)
             )
         ;
 
     // corner bottom right
-    // (should have 5 calls to current concentration)
     row = rowMax-1;
     col = colMax-1;
     concentrations_t1(row,col) = grid.getConcentrations()(row,col) 
         + timestep/(deltaCol*deltaCol)
             * (
-                2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+                calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                    + 2 * grid.getAlphaX()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                * grid.getConcentrations()(row,col-1)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
+                calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaY()(row,col), grid.getAlphaY()(row-1,col))
-                    + 2 * grid.getAlphaY()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaY()(row,col), grid.getAlphaY()(row-1,col))
-                * grid.getConcentrations()(row-1,col)
             )
         ;
 
@@ -238,19 +322,31 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
     return concentrations_t1;
 }
 
-// TODO
-MatrixXd FTCS_closed(Grid grid, Boundary bc, double timestep) {
-    return MatrixXd();
-}
 
 MatrixXd FTCS(Grid grid, Boundary bc, double timestep) {
-    switch (bc.getBoundaryConditionType()) {
+    // inner cells 
-        case BC_TYPE_CONSTANT:
+    // TODO only the boundary cells are different in constant and closed case
-            return FTCS_constant(grid, bc, timestep);
+
-        case BC_TYPE_CLOSED:
+    // if 1D:
-            return FTCS_closed(grid, bc, timestep);
+    //      do inner cells 
-        default:
+    //      do left boundary according to bc type
-            // TODO handle
+    //      do right boundary according to bc type
-            return MatrixXd();
+    // if 2D:
+    //      do inner cells
+    //      do left boundaries according to bc type
+    //      do right boundaries according to bc type
+    //      ...
+
+    if (grid.getDim() == 1) {
+        return FTCS_1D(grid, bc, timestep);
+    } else {
+        return FTCS_2D(grid, bc, timestep);
     }
+
+    // checking the boundary condition type first does not work
+    // if the boundary condition types change dynamically for a grid
+    // meaning:
+    // - boundary condition type needs to be checked for every single boundary cell
+    //      -> this check is last in order
+    // - 
 }

src/Simulation.cpp

@@ -84,7 +84,11 @@ void Simulation::run() {
         }
         printConcentrationsConsole();
         if (csv_output >= CSV_OUTPUT_ON) {
-            printConcentrationsCSV("test", true);
+            bool append = false;
+            if (csv_output == CSV_OUTPUT_VERBOSE) {
+                append = true;
+            }
+            printConcentrationsCSV("test", append);
         }
     } else if (approach == BTCS_APPROACH) {
         for (int i = 0; i < iterations; i++) {

test/CMakeLists.txt

@@ -1,17 +1,20 @@
-include(FetchContent)
+find_library(DOCTEST_LIB doctest)
 
-FetchContent_Declare(
+if(NOT DOCTEST_LIB)
+  include(FetchContent)
+
+  FetchContent_Declare(
     DocTest
     GIT_REPOSITORY https://github.com/doctest/doctest.git
-    GIT_TAG v2.4.9
+    GIT_TAG v2.4.9)
-)
 
-FetchContent_MakeAvailable(DocTest)
+  FetchContent_MakeAvailable(DocTest)
+endif()
 
 add_executable(testTug setup.cpp testBoundaryCondition.cpp testDiffusion.cpp)
 target_link_libraries(testTug doctest tug)
 
-add_custom_target(check
+add_custom_target(
-    COMMAND $<TARGET_FILE:testTug>
+  check
-    DEPENDS testTug
+  COMMAND $<TARGET_FILE:testTug>
-)
+  DEPENDS testTug)

utils/ci.Dockerfile

@@ -0,0 +1,4 @@
+FROM alpine
+    MAINTAINER Max Luebke <mluebke@uni-potsdam.de>
+
+RUN apk add --no-cache build-base openmp cmake git eigen-dev