diff --git a/doc/Implementation.org b/doc/Implementation.org new file mode 100644 index 0000000..d64970b --- /dev/null +++ b/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} diff --git a/include/diffusion/BTCSDiffusion.hpp b/include/diffusion/BTCSDiffusion.hpp index e9eb868..1470bd1 100644 --- a/include/diffusion/BTCSDiffusion.hpp +++ b/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 n_grid_cells Number of grid cells in x direction the domain is + * \param domain_size Size of the domain in x direction. + * \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 n_grid_cells Number of grid cells in y direction the domain is + * \param domain_size Size of the domain in y direction. + * \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 n_grid_cells Number of grid cells in z direction the domain is + * \param domain_size Size of the domain in z direction. + * \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 &c, Eigen::Map &alpha, @@ -131,7 +131,7 @@ private: Eigen::Map &alpha, Eigen::Map &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 &c); diff --git a/scripts/Adi2D_Reference.R b/scripts/Adi2D_Reference.R new file mode 100644 index 0000000..e3c3fd6 --- /dev/null +++ b/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) + diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 5413a63..9d9bd67 100644 --- a/src/BTCSDiffusion.cpp +++ b/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 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), - size, dx, time_step); + fillVectorFromRow(b_vector, c, alpha, bc, d_ortho, size, dx, time_step); // start to solve Eigen::SparseLU, Eigen::COLAMDOrdering> @@ -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 = - calc_t0_c(c.transpose(), alpha.transpose(), bc.transpose(), local_dt, dx); + d_ortho = calc_d_ortho(c.transpose(), alpha.transpose(), bc.transpose(), + 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, - const DMatrixRowMajor &alpha, - const BCMatrixRowMajor &bc, - double time_step, double dx) +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_cols = this->grid_cells[0]; + int n_rows = c.rows(); + int n_cols = c.cols(); - DMatrixRowMajor t0_c(n_rows, n_cols); + DMatrixRowMajor d_ortho(n_rows, n_cols); std::array 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) { - continue; - } - - y_values[0] = getBCFromFlux(tmp_bc, c(0, j), alpha(0, j)); + y_values[0] = (tmp_bc.type == Diffusion::BC_CONSTANT + ? tmp_bc.value + : 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) * - (2 * y_values[0] - 3 * y_values[1] + y_values[2]) / (dx * dx); + d_ortho(0, j) = -sy * (2 * y_values[0] - 3 * y_values[1] + y_values[2]); } // 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) * - (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); + d_ortho(i, j) = -sy * (y_values[0] - 2 * y_values[1] + y_values[2]); } } @@ -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); - - if (tmp_bc.type == Diffusion::BC_CLOSED) { - continue; - } + double sy = (time_step * alpha(end, j)) / (dx * dx); 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) * - (y_values[0] - 3 * y_values[1] + 2 * y_values[2]) / - (dx * dx); + d_ortho(end, j) = -sy * (y_values[0] - 3 * y_values[1] + 2 * y_values[2]); } - 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)); - double value = (c[j] < DOUBLE_MACHINE_EPSILON ? .0 : c[j]); - b_vector[j + 1] = -value - t0_c_j; + b_vector[j + 1] = -c[j] + d_ortho[j]; } // this is not correct currently.We will fix this when we are able to define