Merge branch '8-fix-2d-adi' into 'main'

Resolve "Fix 2D ADI"

Closes #8

See merge request mluebke/diffusion!17

doc/Implementation.org

@@ -0,0 +1,219 @@
+#+title: Implementation of 2D ADI BTCS
+#+LaTeX_CLASS_OPTIONS: [a4paper,10pt]
+#+LATEX_HEADER: \usepackage{fullpage}
+#+LATEX_HEADER: \usepackage{amsmath, systeme}
+#+OPTIONS: toc:nil
+
+* Introduction
+As described in [[./ADI_scheme.org]] the 1D-BTCS diffusion and the 2D-BTCS
+diffusion using the ADI scheme are currently implemented. In the following
+document the implementation is described.
+
+To access columns and rows of input vectors and matrices with ease and to use
+build-in solver algorithms we're using the [[https://eigen.tuxfamily.org/index.php?title=Main_Page][Eigen library]] inside the library.
+However, this not imply to use Eigen as the on-top data structure of the
+application. The API expects pointer to a continuous memory area, where
+concentrations, alpha and boundary conditions are stored. Those memory areas are
+mapped to the according Eigen data structure then.
+
+Solving a linear equation system requires an $A$ matrix and $b$ and $x$ vectors.
+All of them are also implemented using Eigen data structures, whereby the $A$
+matrix is implemented as a [[https://eigen.tuxfamily.org/dox/group__TutorialSparse.html][sparse matrix]].
+
+In the following all indices will start at 0.
+
+** Internal data structures
+
+During instantiating of an object the dimension size is required. According to
+that parameter all other parameters are stored inside the instance of the
+diffusion module and are set with the according getter and setters (see code
+documentation for more details). Stored values are:
+
+- grid dimension
+- time step $\Delta t$
+- number of grid cells $N$ in each direction
+- domain size $x$ or $y$ resp. in each direction
+- spatial discretization in each direction $\Delta x$ or $\Delta y$ resp.
+
+Whenever a value is changed through according setter, the private member
+function =updateInterals()= is called, updating all other dependent values.
+
+** Diffusion in 1D
+
+As described above to call =simulate()= properly, pointers to continuous memory
+containing the field with its concentrations, alphas and boundary conditions are
+required. The boundary condition data type is defined in the
+=BoundayCondition.hpp= header file. There are three assignable types:
+
+- =BC_CONSTANT= - constant/Dirichlet boundary condition with =value= set to a fixed value.
+- =BC_CLOSED= - closed/Neumann boundary condition with =value= set to a gradient of 0
+- =BC_FLUX= - flux/Chauchy boundary with =value= set to user defined gradient
+
+\noindent See the code documentation for more details on that. Defining the
+number of grid cells as $N$, the input values should have the following
+dimensions:
+
+- vector: field =c= $\to N$
+- vector: alpha =alpha= $\to N$
+- vector: boundary condition =bc= $\to N+2$, where =bc[0]= and =bc[N+1]=
+  describing the actual boundary condition on the left/right hand side. All
+  other values can be used to define constant inlet cells.
+
+All input values are mapped to Eigen vectors and passed to the =simulate1D()=
+function. With the only row (since we are just having a 1D vector) we will start
+the simulation by calling =simulate_base()=. There, a =SparseMatrix A= with a
+dimension of $N+1 \times N+1$ and a =VectorXd b= with a dimension of $N+1$ are
+created with, filled with according values and used to solve the LES by the
+built-in =SparseLU= solver of Eigen. The results are written to the =VectorXd x=
+and copied back to the memory area of =c=.
+
+*** Filling of sparse matrix =A=
+
+To fill the =A= matrix various information about the field are needed. This includes:
+
+- =alpha= for each cell of the field
+- =bc= for each cell of the field and boundary conditions
+- =size= of the field
+- =dx= as the $\Delta x$
+- =time_step= as the current time step size $\Delta t$
+
+We are using indices $j$ from $0 \to N+1$ to address the correct field of the
+boundary condition and $k$ from $0 \to N-1$ to address the according field of
+=alpha=
+
+First of all the left and right boundary condition are extracted from =bc=
+evaluated. Then we will set the left boundary condition to the linear equation
+system, which is the left boundary condition, to 1 ($A(0,0) = 1$). Now solving
+for the leftmost inlet is required. Each $sx = (\alpha_k \cdot \Delta t)/ \Delta
+x^2$ is stored in the same named variable for each index pointing to the current
+cell of the input/output field.
+
+- $A(1,0) = 2 \cdot sx$
+- $A(1,1) = -1 -3 \cdot sx$
+- $A(1,2) = sx$
+
+After that we will iterate until $k = N-2$ or $j = N-1$ is reached respectively
+and setting =A= with =alpha[k]= and =bc[j]= to the following values:
+
+- $A(j,j-1) = sx$
+- $A(j,j) = -1 -2 \cdot sx$
+- $A(j,j+1) = sx$
+
+Now the equation for the rightmost inlet cell is applied:
+
+- $A(N,N-1) = sx$
+- $A(N,N) = -1 -3 \cdot sx$
+- $A(N,N+1) = 2 \cdot sx$
+
+The right boundary condition is also set to $A(N+1,N+1) = 1$.
+
+A special case occurs if one cell of the input/output field is set to a constant
+value by setting the according =bc= field to =BC_CONSTANT=. Then, just the according
+element of the matrix is set to 1. For example cell $k=4$ is marked as
+constant will lead to $j = 5$ and $A(j,j) = 1$.
+
+*** Filling of vector =b=
+
+The right hand side of the equation system requires the following input values:
+
+- =c= as the input field
+- =alpha= for each cell of the grid
+- =bc= for each cell of the field and boundary conditions
+- =d_ortho= as the explicit parameter of left hand side of the equation for each cell of the field
+  (only interesting if using 2D ADI, otherwise this parameter is set to 0 for
+  all cells)
+- =size= of the field
+- =dx= as the $\Delta x$
+- =time_step= as the current time step size $\Delta t$
+
+Now we are able to iterating through the indices of $j$ from $0 \to N-1$ and
+filling up the =b= vector by applying the following rule:
+
+\begin{equation}\displaystyle
+b_{j+1} = \begin{cases}
+bc_{j+1}.\text{value} & bc_{j+1}.\text{type} = \texttt{BC\_CONSTANT} \\
+-\texttt{c}_j - \texttt{d\_ortho}_j & \text{else}
+\end{cases}
+\end{equation}
+
+\noindent With
+
+\begin{equation}\displaystyle
+\texttt{d\_ortho}_j = \Delta t \cdot \alpha_j \cdot (\frac{\bot d_j}{\Delta x^2})
+\end{equation}
+
+Both boundary conditions on the left and right side are assigned by either using
+a constant value or a derived value from a given gradient, depending on the
+assigned boundary condition type. Until now, only gradients of 0 (closed
+boundary condition) are supported. This leads to:
+
+\begin{equation}\displaystyle
+b_{0} = \begin{cases}
+bc_0.\text{value} & bc_0.\text{type} = \texttt{BC\_CONSTANT} \\
+c_0 & \text{else}
+\end{cases}
+\end{equation}
+
+\begin{equation}\displaystyle
+b_{N+1} = \begin{cases}
+bc_{N+1}.\text{value} & bc_{N+1}.\text{type} = \texttt{BC\_CONSTANT} \\
+c_{N-1} & \text{else}
+\end{cases}
+\end{equation}
+** Diffusion in 2D
+
+Simulating a 2D diffusion works in overall just like the 1D diffusion. Input
+parameters are should now be matrices or more precisly arrays, which can be
+interpreted as matrices, no matter if stored column or row wise. Defining the
+number of grid cells now as $N times M$, the input parameters should have the
+following dimensions:
+
+- vector: field =c= $\to N \times M$
+- vector: alpha =alpha= $\to N \times M$
+- vector: boundary condition =bc= $\to N+2 \times M+2$
+
+Those input pointers will now be interpreted as Eigen matrices for further use.
+In =simulate2D()= the matrix =d_ortho= is calculated in y-direction. Then, for
+each row of matrix =c= =simulate_base()= is called. The according row of
+=alpha=, =bc= and =d_ortho= are passed too and the time step is divided by two.
+Now the 2D problem is divided into an 1D problem. After each row the new
+resulting field is written back into the original memory area and serves as the
+input of the second simulation in the y-direction.
+
+For this, =d_ortho= will be recalculated, now in x-direction. Now, for every
+column of =c=, the correspondending column =alpha=, =bc= and =d_ortho= is passed
+to =simulate_base()= with also half the time step.
+
+Finally, the result of the simulation of each column get written back to =c= and
+the 2D simulation is finished.
+
+*** Calculation of =d_ortho=
+
+The function =cal_d_ortho= returns a matrix with the size of $N \times M$ with
+each cell holding the $\bot d$ value of a cell from =c= at the same index. This
+is done by first iterating over the first row from $0 \to N-1$. There, we apply
+the following rule:
+
+\begin{equation}\displaystyle
+\bot d_{0,j} = \begin{cases}
+\Delta t \cdot \alpha_{0,j} \cdot \left(\frac{2bc_{0,j+1}.\text{value} - 3 \cdot c_{0,j} + c_{1,j}}{\Delta x^2}\right) & bc_{0,j+1} = \texttt{BC\_CONSTANT} \\
+\Delta t \cdot \alpha_{0,j} \cdot \left(\frac{c_{0,j} - 2 \cdot c_{0,j} + c_{1,j}}{\Delta x^2}\right) & \text{else}
+\end{cases}
+\end{equation}
+
+\noindent Then we succeeding to the inner cells with $i = 1,\dots,M-2$ and $j =
+0,\dots,N-1$:
+
+\begin{equation}\displaystyle
+\bot d_{i,j} = \Delta t \cdot \alpha_{i,j} \cdot \left(\frac{c_{i-1,j} - 2 \cdot c_{i,j} + c_{i+1,j}}{\Delta x^2}\right)
+\end{equation}
+
+\noindent And lastly calculating the last row by applying the following rule for
+each $j = 0,\dots,N-1$:
+
+\begin{equation}\displaystyle
+\bot d_{M-1,j} = \begin{cases}
+\Delta t \cdot \alpha_{M-1,j} \cdot \left(c_{M-2,j} -3 \cdot c_{M-1,j} + 2\frac{bc_{M+1,j+1}.\text{value}}{\Delta x^2}\right) & bc_{M+1,j+1} = \texttt{BC\_CONSTANT} \\
+\Delta t \cdot \alpha_{M-1,j} \cdot \left(\frac{c_{M-2,j} - 2 \cdot c_{M-1,j} + c_{M-1,j}}{\Delta x^2}\right) & \text{else}
+\end{cases}
+\end{equation}

include/diffusion/BTCSDiffusion.hpp

@@ -24,7 +24,7 @@ public:
   /*!
    * Creates a diffusion module.
    *
-   * @param dim Number of dimensions. Should not be greater than 3 and not less
+   * \param dim Number of dimensions. Should not be greater than 3 and not less
    * than 1.
    */
   BTCSDiffusion(unsigned int dim);
@@ -32,8 +32,8 @@ public:
   /*!
    * Define the grid in x direction.
    *
-   * @param domain_size Size of the domain in x direction.
+   * \param domain_size Size of the domain in x direction.
-   * @param n_grid_cells Number of grid cells in x direction the domain is
+   * \param n_grid_cells Number of grid cells in x direction the domain is
    * divided to.
    */
   void setXDimensions(double domain_size, unsigned int n_grid_cells);
@@ -43,8 +43,8 @@ public:
    *
    * Throws an error if the module wasn't initialized at least as a 2D model.
    *
-   * @param domain_size Size of the domain in y direction.
+   * \param domain_size Size of the domain in y direction.
-   * @param n_grid_cells Number of grid cells in y direction the domain is
+   * \param n_grid_cells Number of grid cells in y direction the domain is
    * divided to.
    */
   void setYDimensions(double domain_size, unsigned int n_grid_cells);
@@ -54,8 +54,8 @@ public:
    *
    * Throws an error if the module wasn't initialized at least as a 3D model.
    *
-   * @param domain_size Size of the domain in z direction.
+   * \param domain_size Size of the domain in z direction.
-   * @param n_grid_cells Number of grid cells in z direction the domain is
+   * \param n_grid_cells Number of grid cells in z direction the domain is
    * divided to.
    */
   void setZDimensions(double domain_size, unsigned int n_grid_cells);
@@ -89,13 +89,13 @@ public:
   /*!
    * With given ghost zones simulate diffusion. Only 1D allowed at this moment.
    *
-   * @param c Pointer to continious memory describing the current concentration
+   * \param c Pointer to continious memory describing the current concentration
    * state of each grid cell.
-   * @param alpha Pointer to memory area of diffusion coefficients for each grid
+   * \param alpha Pointer to memory area of diffusion coefficients for each grid
    * element.
-   * @param bc Pointer to memory setting boundary conidition of each grid cell.
+   * \param bc Pointer to memory setting boundary conidition of each grid cell.
    *
-   * @return Time in seconds [s] used to simulate one iteration.
+   * \return Time in seconds [s] used to simulate one iteration.
    */
   auto simulate(double *c, double *alpha, Diffusion::boundary_condition *bc)
       -> double;
@@ -103,7 +103,7 @@ public:
   /*!
    * Set the timestep of the simulation
    *
-   * @param time_step Time step (in seconds ???)
+   * \param time_step Time step (in seconds ???)
    */
   void setTimestep(double time_step);
 
@@ -121,7 +121,7 @@ private:
 
   void simulate_base(DVectorRowMajor &c, const BCVectorRowMajor &bc,
                      const DVectorRowMajor &alpha, double dx, double time_step,
-                     int size, const DVectorRowMajor &t0_c);
+                     int size, const DVectorRowMajor &d_ortho);
 
   void simulate1D(Eigen::Map<DVectorRowMajor> &c,
                   Eigen::Map<const DVectorRowMajor> &alpha,
@@ -131,7 +131,7 @@ private:
                   Eigen::Map<const DMatrixRowMajor> &alpha,
                   Eigen::Map<const BCMatrixRowMajor> &bc);
 
-  auto calc_t0_c(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
+  auto calc_d_ortho(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
                  const BCMatrixRowMajor &bc, double time_step, double dx)
       -> DMatrixRowMajor;
 
@@ -144,7 +144,7 @@ private:
                                 const DVectorRowMajor &c,
                                 const DVectorRowMajor &alpha,
                                 const BCVectorRowMajor &bc,
-                                const DVectorRowMajor &t0_c, int size,
+                                const DVectorRowMajor &d_ortho, int size,
                                 double dx, double time_step);
   void simulate3D(std::vector<double> &c);
 

scripts/Adi2D_Reference.R

@@ -0,0 +1,124 @@
+
+## Brutal implementation of 2D ADI scheme
+## Square NxN grid with dx=dy=1
+ADI <- function(n, dt, iter, alpha) {
+
+    nx <- ny <- n
+    dx <- dy <- 1
+    field <- matrix(0, nx, ny)
+    ## find out the center of the grid to apply conc=1
+    cen <- ceiling(n/2)
+    field[cen, cen] <- 1
+
+    ## prepare containers for computations and outputs
+    tmpX <- tmpY <- res <- field
+    out <- vector(mode="list", length=iter)
+    
+    for (it in seq(1, iter)) {
+        for (i in seq(1, ny))
+            tmpX[i,] <- SweepByRow(i, res, dt=dt, alpha=alpha)
+
+        resY <- t(tmpX)
+        for (i in seq(1, nx))
+            tmpY[i,] <- SweepByRow(i, resY, dt=dt, alpha=alpha)
+
+        res <- t(tmpY)
+        out[[it]] <- res }
+
+    return(out)
+}
+
+
+## Workhorse function to fill A, B and solve for a given *row* of the
+## grid matrix
+SweepByRow <- function(i, field, dt, alpha) {
+    dx <- 1 ## fixed in our test
+    A <- matrix(0, nrow(field), ncol(field))
+    Sx <- Sy <- alpha*dt/2/dx/dx
+
+    ## diagonal of A at once
+    diag(A) <- -1-2*Sx
+
+    ## adjacent diagonals "Sx"
+    for (ii in seq(1, nrow(field)-1)){
+        A[ii+1, ii] <- Sx
+        A[ii, ii+1] <- Sx
+    }
+
+    B <- numeric(ncol(field))
+
+    ## We now distinguish the top and bottom rows
+    if (i == 1){
+        ## top boundary, "i-1" doesn't exist or is at a ghost
+        ## node/cell boundary (TODO)
+        for (ii in seq_along(B))
+            B[ii] <- (-1 +2*Sy)*field[i,ii] - Sy*field[i+1,ii]
+    } else if (i == nrow(field)){ 
+        ## bottom boundary, "i+1" doesn't exist or is at a ghost
+        ## node/cell boundary (TODO)
+        for (ii in seq_along(B))
+            B[ii] <- -Sy*field[i-1, ii] + (-1 +2*Sy)*field[i,ii]
+        
+    } else {
+        ## inner grid row, full expression
+        for (ii in seq_along(B))
+            B[ii] <- -Sy*field[i-1, ii] + (-1 +2*Sy)*field[i,ii] - Sy*field[i+1,ii]
+    }
+    
+    x <- solve(A, B)
+    x
+}
+
+
+
+DoRef <- function(n, alpha, dt, iter) {
+    require(ReacTran)
+    require(deSolve)
+    
+    N     <- n          # number of grid cells
+    XX    <- n          # total size
+    dy    <- dx <- XX/N  # grid size
+    Dy    <- Dx <- alpha  # diffusion coeff, X- and Y-direction
+    
+    ## The model equations
+    
+    Diff2D <- function (t, y, parms)  {
+        CONC  <- matrix(nrow = N, ncol = N, y)
+        dCONC <- tran.2D(CONC, D.x = Dx, D.y = Dy, dx = dx, dy = dy)$dC
+        return (list(dCONC))
+    }
+    
+    ## initial condition: 0 everywhere, except in central point
+    y <- matrix(nrow = N, ncol = N, data = 0)
+    cen <- ceiling(N/2)
+    y[cen, cen] <- 1  ## initial concentration in the central point...
+    
+
+    ## solve for  time units
+    times <- seq(0,iter)*dt
+    
+    out <- ode.2D (y = y, func = Diff2D, t = times, parms = NULL,
+                   dim = c(N,N), lrw=155412)
+
+    ref <- matrix(out[length(times),-1], N, N)
+    return(ref)
+}
+
+## test number 1
+adi1 <- ADI(n=25, dt=100, iter=50, alpha=1E-3)
+ref1 <- DoRef(n=25, alpha=1E-3, dt=100, iter=50)
+
+plot(adi1[[length(adi1)]], ref1, log="xy", xlab="ADI", ylab="ode.2D (reference)",
+     las=1, xlim=c(1E-15, 1), ylim=c(1E-15, 1))
+abline(0,1)
+
+sapply(adi1, sum)
+
+## test number 2
+adi2 <- ADI(n=51, dt=10, iter=200, alpha=1E-3)
+ref2 <- DoRef(n=51, alpha=1E-3, dt=10, iter=200)
+
+plot(adi2[[length(adi2)]], ref2, log="xy", xlab="ADI", ylab="ode.2D (reference)",
+     las=1, xlim=c(1E-15, 1), ylim=c(1E-15, 1))
+abline(0,1)
+

src/BTCSDiffusion.cpp

@@ -91,7 +91,7 @@ void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c,
                                              const DVectorRowMajor &alpha,
                                              double dx, double time_step,
                                              int size,
-                                             const DVectorRowMajor &t0_c) {
+                                             const DVectorRowMajor &d_ortho) {
 
   Eigen::SparseMatrix<double> A_matrix;
   Eigen::VectorXd b_vector;
@@ -104,8 +104,7 @@ void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c,
   x_vector.resize(size + 2);
 
   fillMatrixFromRow(A_matrix, alpha, bc, size, dx, time_step);
-  fillVectorFromRow(b_vector, c, alpha, bc, Eigen::VectorXd::Constant(size, 0),
+  fillVectorFromRow(b_vector, c, alpha, bc, d_ortho, size, dx, time_step);
-                    size, dx, time_step);
 
   // start to solve
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
@@ -145,70 +144,68 @@ void Diffusion::BTCSDiffusion::simulate2D(
 
   double local_dt = this->time_step / BTCS_2D_DT_SIZE;
 
-  DMatrixRowMajor t0_c = calc_t0_c(c, alpha, bc, local_dt, dx);
+  DMatrixRowMajor d_ortho = calc_d_ortho(c, alpha, bc, local_dt, dx);
 
 #pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < n_rows; i++) {
     DVectorRowMajor input_field = c.row(i);
     simulate_base(input_field, bc.row(i + 1), alpha.row(i), dx, local_dt,
-                  n_cols, t0_c.row(i));
+                  n_cols, d_ortho.row(i));
     c.row(i) << input_field;
   }
 
   dx = this->deltas[1];
 
-  t0_c =
+  d_ortho = calc_d_ortho(c.transpose(), alpha.transpose(), bc.transpose(),
-      calc_t0_c(c.transpose(), alpha.transpose(), bc.transpose(), local_dt, dx);
+                         local_dt, dx);
 
 #pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < n_cols; i++) {
     DVectorRowMajor input_field = c.col(i);
     simulate_base(input_field, bc.col(i + 1), alpha.col(i), dx, local_dt,
-                  n_rows, t0_c.row(i));
+                  n_rows, d_ortho.row(i));
     c.col(i) << input_field.transpose();
   }
 }
 
-auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
+auto Diffusion::BTCSDiffusion::calc_d_ortho(const DMatrixRowMajor &c,
                                             const DMatrixRowMajor &alpha,
                                             const BCMatrixRowMajor &bc,
                                             double time_step, double dx)
     -> DMatrixRowMajor {
 
-  int n_rows = this->grid_cells[1];
+  int n_rows = c.rows();
-  int n_cols = this->grid_cells[0];
+  int n_cols = c.cols();
 
-  DMatrixRowMajor t0_c(n_rows, n_cols);
+  DMatrixRowMajor d_ortho(n_rows, n_cols);
 
   std::array<double, 3> y_values{};
 
   // first, iterate over first row
   for (int j = 0; j < n_cols; j++) {
     boundary_condition tmp_bc = bc(0, j + 1);
+    double sy = (time_step * alpha(0, j)) / (dx * dx);
 
-    if (tmp_bc.type == Diffusion::BC_CLOSED) {
+    y_values[0] = (tmp_bc.type == Diffusion::BC_CONSTANT
-      continue;
+                       ? tmp_bc.value
-    }
+                       : getBCFromFlux(tmp_bc, c(0, j), alpha(0, j)));
-
-    y_values[0] = getBCFromFlux(tmp_bc, c(0, j), alpha(0, j));
     y_values[1] = c(0, j);
     y_values[2] = c(1, j);
 
-    t0_c(0, j) = time_step * alpha(0, j) *
+    d_ortho(0, j) = -sy * (2 * y_values[0] - 3 * y_values[1] + y_values[2]);
-                 (2 * y_values[0] - 3 * y_values[1] + y_values[2]) / (dx * dx);
   }
 
 // then iterate over inlet
 #pragma omp parallel for private(y_values) schedule(dynamic)
   for (int i = 1; i < n_rows - 1; i++) {
     for (int j = 0; j < n_cols; j++) {
+      double sy = (time_step * alpha(i, j)) / (dx * dx);
 
       y_values[0] = c(i - 1, j);
       y_values[1] = c(i, j);
       y_values[2] = c(i + 1, j);
 
-      t0_c(i, j) = time_step * alpha(i, j) *
+      d_ortho(i, j) = -sy * (y_values[0] - 2 * y_values[1] + y_values[2]);
-                   (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx);
     }
   }
 
@@ -217,21 +214,18 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c,
   // and finally over last row
   for (int j = 0; j < n_cols; j++) {
     boundary_condition tmp_bc = bc(end + 1, j + 1);
-
+    double sy = (time_step * alpha(end, j)) / (dx * dx);
-    if (tmp_bc.type == Diffusion::BC_CLOSED) {
-      continue;
-    }
 
     y_values[0] = c(end - 1, j);
     y_values[1] = c(end, j);
-    y_values[2] = getBCFromFlux(tmp_bc, c(end, j), alpha(end, j));
+    y_values[2] = (tmp_bc.type == Diffusion::BC_CONSTANT
+                       ? tmp_bc.value
+                       : getBCFromFlux(tmp_bc, c(end, j), alpha(end, j)));
 
-    t0_c(end, j) = time_step * alpha(end, j) *
+    d_ortho(end, j) = -sy * (y_values[0] - 3 * y_values[1] + 2 * y_values[2]);
-                   (y_values[0] - 3 * y_values[1] + 2 * y_values[2]) /
-                   (dx * dx);
   }
 
-  return t0_c;
+  return d_ortho;
 }
 
 void Diffusion::BTCSDiffusion::fillMatrixFromRow(
@@ -285,7 +279,7 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow(
 void Diffusion::BTCSDiffusion::fillVectorFromRow(
     Eigen::VectorXd &b_vector, const DVectorRowMajor &c,
     const DVectorRowMajor &alpha, const BCVectorRowMajor &bc,
-    const DVectorRowMajor &t0_c, int size, double dx, double time_step) {
+    const DVectorRowMajor &d_ortho, int size, double dx, double time_step) {
 
   Diffusion::boundary_condition left = bc[0];
   Diffusion::boundary_condition right = bc[size + 1];
@@ -303,9 +297,7 @@ void Diffusion::BTCSDiffusion::fillVectorFromRow(
       continue;
     }
 
-    double t0_c_j = time_step * alpha[j] * (t0_c[j] / (dx * dx));
+    b_vector[j + 1] = -c[j] + d_ortho[j];
-    double value = (c[j] < DOUBLE_MACHINE_EPSILON ? .0 : c[j]);
-    b_vector[j + 1] = -value - t0_c_j;
   }
 
   // this is not correct currently.We will fix this when we are able to define