diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 2e36801..e5f57a5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -19,7 +19,7 @@ build: expire_in: 100s script: - mkdir build && cd build - - cmake .. + - cmake -DCMAKE_BUILD_TYPE=Debug .. - make run_1D: @@ -41,7 +41,7 @@ lint: script: - mkdir lint && cd lint - cmake -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_CXX_CLANG_TIDY="clang-tidy;-checks=cppcoreguidelines-*,clang-analyzer-*,performance-*,readability-*, modernize-*" .. - - make + - make diffusion memcheck_1D: stage: dynamic_analyze diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt index c92d35b..62b0bfd 100644 --- a/app/CMakeLists.txt +++ b/app/CMakeLists.txt @@ -6,3 +6,6 @@ target_link_libraries(2D PUBLIC diffusion) add_executable(Comp2D main_2D_mdl.cpp) target_link_libraries(Comp2D PUBLIC diffusion) + +add_executable(Const2D main_2D_const.cpp) +target_link_libraries(Const2D PUBLIC diffusion) diff --git a/app/main_1D.cpp b/app/main_1D.cpp index b2423f2..2fb7c21 100644 --- a/app/main_1D.cpp +++ b/app/main_1D.cpp @@ -20,7 +20,7 @@ int main(int argc, char *argv[]) { // create input + diffusion coefficients for each grid cell std::vector alpha(n, 1 * pow(10, -1)); std::vector field(n, 1 * std::pow(10, -6)); - std::vector bc(n, {0,0}); + std::vector bc(n + 2, {0, 0}); // create instance of diffusion module BTCSDiffusion diffu(dim); @@ -28,7 +28,7 @@ int main(int argc, char *argv[]) { diffu.setXDimensions(1, n); // set the boundary condition for the left ghost cell to dirichlet - bc[0] = {Diffusion::BC_CONSTANT, 5*std::pow(10,-6)}; + bc[0] = {Diffusion::BC_CONSTANT, 5 * std::pow(10, -6)}; // diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT, // 5. * std::pow(10, -6)); @@ -44,7 +44,7 @@ int main(int argc, char *argv[]) { cout << "Iteration: " << i << "\n\n"; - for (auto & cell : field) { + for (auto &cell : field) { cout << cell << "\n"; } diff --git a/app/main_2D.cpp b/app/main_2D.cpp index 6a35bba..502ad0e 100644 --- a/app/main_2D.cpp +++ b/app/main_2D.cpp @@ -1,7 +1,7 @@ +#include // for copy, max +#include #include #include -#include // for copy, max -#include #include #include // for std #include // for vector @@ -20,7 +20,7 @@ int main(int argc, char *argv[]) { // create input + diffusion coefficients for each grid cell std::vector alpha(n * m, 1 * pow(10, -1)); std::vector field(n * m, 1 * std::pow(10, -6)); - std::vector bc(n*m, {0,0}); + std::vector bc((n + 2) * (m + 2), {0, 0}); // create instance of diffusion module BTCSDiffusion diffu(dim); @@ -28,8 +28,8 @@ int main(int argc, char *argv[]) { diffu.setXDimensions(1, n); diffu.setYDimensions(1, m); - //set inlet to higher concentration without setting bc - field[12] = 5*std::pow(10,-3); + // set inlet to higher concentration without setting bc + field[12] = 5 * std::pow(10, -3); // set timestep for simulation to 1 second diffu.setTimestep(1.); diff --git a/app/main_2D_const.cpp b/app/main_2D_const.cpp new file mode 100644 index 0000000..e70f8ad --- /dev/null +++ b/app/main_2D_const.cpp @@ -0,0 +1,57 @@ +#include // for copy, max +#include +#include +#include +#include +#include // for std +#include // for vector +using namespace std; + +using namespace Diffusion; + +int main(int argc, char *argv[]) { + + // dimension of grid + int dim = 2; + + int n = 5; + int m = 5; + + // create input + diffusion coefficients for each grid cell + std::vector alpha(n * m, 1 * pow(10, -1)); + std::vector field(n * m, 1 * std::pow(10, -6)); + std::vector bc((n + 2) * (m + 2), {0, 0}); + + // create instance of diffusion module + BTCSDiffusion diffu(dim); + + diffu.setXDimensions(1, n); + diffu.setYDimensions(1, m); + + for (int i = 1; i <= n; i++) { + bc[(n + 2) * i] = {Diffusion::BC_CONSTANT, 5 * std::pow(10, -6)}; + } + // set timestep for simulation to 1 second + diffu.setTimestep(1.); + + cout << setprecision(12); + + for (int t = 0; t < 10; t++) { + diffu.simulate(field.data(), alpha.data(), bc.data()); + + cout << "Iteration: " << t << "\n\n"; + + // iterate through rows + for (int i = 0; i < m; i++) { + // iterate through columns + for (int j = 0; j < n; j++) { + cout << left << std::setw(20) << field[i * n + j]; + } + cout << "\n"; + } + + cout << "\n" << endl; + } + + return 0; +} diff --git a/app/main_2D_mdl.cpp b/app/main_2D_mdl.cpp index c040cfe..f4eab0c 100644 --- a/app/main_2D_mdl.cpp +++ b/app/main_2D_mdl.cpp @@ -1,7 +1,7 @@ +#include // for copy, max +#include #include #include -#include // for copy, max -#include #include #include // for std #include // for vector @@ -20,7 +20,7 @@ int main(int argc, char *argv[]) { // create input + diffusion coefficients for each grid cell std::vector alpha(n * m, 1 * pow(10, -1)); std::vector field(n * m, 0.); - std::vector bc(n*m, {0,0}); + std::vector bc((n + 2) * (m + 2), {0, 0}); field[125500] = 1; @@ -31,13 +31,13 @@ int main(int argc, char *argv[]) { diffu.setYDimensions(1., m); // set the boundary condition for the left ghost cell to dirichlet - //diffu.setBoundaryCondition(250, 250, BTCSDiffusion::BC_CONSTANT, 1); + // diffu.setBoundaryCondition(250, 250, BTCSDiffusion::BC_CONSTANT, 1); // for (int d=0; d<5;d++){ // diffu.setBoundaryCondition(d, 0, BC_CONSTANT, .1); // } // diffu.setBoundaryCondition(1, 1, BTCSDiffusion::BC_CONSTANT, .1); // diffu.setBoundaryCondition(1, 1, BTCSDiffusion::BC_CONSTANT, .1); - + // set timestep for simulation to 1 second diffu.setTimestep(1.); @@ -45,9 +45,9 @@ int main(int argc, char *argv[]) { // First we output the initial state cout << 0; - - for (int i=0; i < m*n; i++) { - cout << "," << field[i]; + + for (int i = 0; i < m * n; i++) { + cout << "," << field[i]; } cout << endl; @@ -58,9 +58,9 @@ int main(int argc, char *argv[]) { cerr << "time elapsed: " << time << " seconds" << endl; cout << t; - - for (int i=0; i < m*n; i++) { - cout << "," << field[i]; + + for (int i = 0; i < m * n; i++) { + cout << "," << field[i]; } cout << endl; } diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 6a86caf..50544fe 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -1,9 +1,314 @@ -#+TITLE: Adi 2D Scheme +#+TITLE: Numerical solution of diffusion equation in 2D with ADI Scheme #+LaTeX_CLASS_OPTIONS: [a4paper,10pt] #+LATEX_HEADER: \usepackage{fullpage} +#+LATEX_HEADER: \usepackage{amsmath, systeme} #+OPTIONS: toc:nil -* Input + +* Diffusion in 1D + +** Finite differences with nodes as cells' centres + +The 1D diffusion equation is: + +#+NAME: eqn:1 +\begin{align} +\frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha \frac{\partial C }{\partial x} \right) \nonumber \\ + & = \alpha \frac{\partial^2 C}{\partial x^2} +\end{align} + +We aim at numerically solving [[eqn:1]] on a spatial grid such as: + +[[./grid_pqc.pdf]] + +The left boundary is defined on $x=0$ while the center of the first +cell - which are the points constituting the finite difference nodes - +is in $x=dx/2$, with $dx=L/n$. + + +** The explicit FTCS scheme (as in PHREEQC) + +We start by discretizing [[eqn:1]] following an explicit Euler scheme and +specifically a Forward Time, Centered Space finite difference. + +For each cell index $i \in 1, \dots, n-1$ and assuming constant +$\alpha$, we can write: + +#+NAME: eqn:2 +\begin{equation}\displaystyle + \frac{C_i^{t+1} -C_i^{t}}{\Delta t} = \alpha\frac{\frac{C^t_{i+1}-C^t_{i}}{\Delta x}-\frac{C^t_{i}-C^t_{i-1}}{\Delta x}}{\Delta x} +\end{equation} + +In practice, we evaluate the first derivatives of $C$ w.r.t. $x$ on +the boundaries of each cell (i.e., $(C_{i+1}-C_i)/\Delta x$ on the +right boundary of the i-th cell and $(C_{i}-C_{i-1})/\Delta x$ on its +left cell boundary) and then repeat the differentiation to get the +second derivative of $C$ on the the cell centre $i$. + +This discretization works for all internal cells, but not for the +domain boundaries ($i=0$ and $i=n$). To properly treat them, we need +to account for the discrepancy in the discretization. + +For the first (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: + +#+NAME: eqn:3 +\begin{equation}\displaystyle +\frac{C_0^{t+1} -C_0^{t}}{\Delta t} = \alpha\frac{\frac{C^t_{1}-C^t_{0}}{\Delta x}- +\frac{C^t_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +#+NAME: eqn:4 +\begin{align}\displaystyle +C_0^{t+1} & = C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^t_{1}-C^t_{0}- 2 C^t_{0}+2l \right) \nonumber \\ + & = C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^t_{1}- 3 C^t_{0} +2l \right) +\end{align} + + +In case of constant right boundary, the finite difference of point +$C_n$ - calling $r$ the right boundary value - is: + +#+NAME: eqn:5 +\begin{equation}\displaystyle +\frac{C_n^{t+1} -C_n^t}{\Delta t} = \alpha\frac{\frac{r - C^t_{n}}{\frac{\Delta x}{2}}- +\frac{C^t_{n}-C^t_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +Which, developed, gives +#+NAME: eqn:6 +\begin{align}\displaystyle +C_n^{t+1} & = C_n^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^t_{n} -C^t_{n} + C^t_{n-1} \right) \nonumber \\ + & = C_n^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^t_{n} + C^t_{n-1} \right) +\end{align} + +If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:5]] +becomes zero and we are left with: + + +#+NAME: eqn:7 +\begin{equation}\displaystyle +C_n^{t+1} = C_n^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^t_{n-1} - C^t_n) +\end{equation} + + + +A similar treatment can be applied to the BTCS implicit scheme. + +** Implicit BTCS scheme + +First, we define the Backward Time difference: + +\begin{equation} + \frac{\partial C^{t+1} }{\partial t} = \frac{C^{t+1}_i - C^{t}_i}{\Delta t} +\end{equation} + +Second the spatial derivative approximation, evaluated at time level $t+1$: + +#+NAME: eqn:secondderivative +\begin{equation} + \frac{\partial^2 C^{t+1} }{\partial x^2} = \frac{\frac{C^{t+1}_{i+1}-C^{t+1}_{i}}{\Delta x}-\frac{C^{t+1}_{i}-C^{t+1}_{i-1}}{\Delta x}}{\Delta x} +\end{equation} + +Taking the 1D diffusion equation from [[eqn:1]] and substituting each term by the +equations given above leads to the following equation: + + +# \begin{equation}\displaystyle +# \frac{C_i^{j+1} -C_i^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{i+1}-C^{j+1}_{i}}{\Delta x}-\frac{C^{j+1}_{i}-C^{j+1}_{i-1}}{\Delta x}}{\Delta x} +# \end{equation} + +# Since we are not able to solve this system w.r.t unknown values in $C^{j-1}$ we +# are shifting each j by 1 to $j \to (j+1)$ and $(j-1) \to j$ which leads to: + +#+NAME: eqn:1DBTCS +\begin{align}\displaystyle +\frac{C_i^{t+1} - C_i^{t}}{\Delta t} & = \alpha\frac{\frac{C^{t+1}_{i+1}-C^{t+1}_{i}}{\Delta x}-\frac{C^{t+1}_{i}-C^{t+1}_{i-1}}{\Delta x}}{\Delta x} \nonumber \\ + & = \alpha\frac{C^{t+1}_{i-1} - 2C^{t+1}_{i} + C^{t+1}_{i+1}}{\Delta x^2} +\end{align} + +This only applies to inlet cells with no ghost node as neighbor. For the left +cell with its center at $\frac{dx}{2}$ and the constant concentration on the +left ghost node called $l$ the equation goes as followed: + +\begin{equation}\displaystyle +\frac{C_0^{t+1} -C_0^{t}}{\Delta t} = \alpha\frac{\frac{C^{t+1}_{1}-C^{t+1}_{0}}{\Delta x}- +\frac{C^{t+1}_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +\begin{align}\displaystyle +C_0^{t+1} & = C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{t+1}_{1}-C^{t+1}_{0}- 2 C^{t+1}_{0}+2l \right) \nonumber \\ + & = C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{t+1}_{1}- 3 C^{t+1}_{0} +2l \right) +\end{align} + +Now we define variable $s_x$ as followed: + +\begin{equation} + s_x = \frac{\alpha \cdot \Delta t}{\Delta x^2} +\end{equation} + +Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: + +\begin{equation}\displaystyle + -C^t_0 = (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{t+1}_0 + s_x \cdot C^{t+1}_1 +\end{equation} + +The right boundary follows the same scheme. We now want to show the equation for the rightmost inlet cell $C_n$ with right boundary value $r$: + +\begin{equation}\displaystyle +\frac{C_n^{t+1} -C_n^{t}}{\Delta t} = \alpha\frac{\frac{r-C^{t+1}_{n}}{\frac{\Delta x}{2}}- +\frac{C^{t+1}_{n}-C^{t+1}_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +\begin{align}\displaystyle +C_n^{t+1} & = C_n^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 2C^{t+1}_{n} - C^{t+1}_{n} + C^{t+1}_{n-1} \right) \nonumber \\ + & = C_0^{t} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 3C^{t+1}_{n} + C^{t+1}_{n-1} \right) +\end{align} + +Now rearrange terms and substituting with $s_x$ leads to: + +\begin{equation}\displaystyle + -C^t_n = s_x \cdot C^{t+1}_{n-1} + (-1 - 3s_x) \cdot C^{t+1}_n + (2s_x) \cdot r +\end{equation} + +*TODO* +- Tridiagonal matrix filling + +#+LATEX: \clearpage + +* Diffusion in 2D: the Alternating Direction Implicit scheme + + +In 2D, the diffusion equation (in absence of source terms) and +assuming homogeneous but anisotropic diffusion coefficient +$\alpha=(\alpha_x,\alpha_y)$ becomes: + +#+NAME: eqn:2d +\begin{equation} +\displaystyle \frac{\partial C}{\partial t} = \alpha_x \frac{\partial^2 C}{\partial x^2} + \alpha_y\frac{\partial^2 C}{\partial y^2} +\end{equation} + +** 2D ADI using BTCS scheme + +The Alternating Direction Implicit method consists in splitting the +integration of eq. [[eqn:2d]] in two half-steps, each of which represents +implicitly the derivatives in one direction, and explicitly in the +other. Therefore we make use of second derivative operator defined in +equation [[eqn:secondderivative]] in both $x$ and $y$ direction for half +the time step $\Delta t$. + +Denoting $i$ the grid cell index along $x$ direction, $j$ the index in +$y$ direction and $t$ as the time level, the spatially centered second +derivatives can be written as: + +\begin{align}\displaystyle +\frac{\partial^2 C^t_{i,j}}{\partial x^2} &= \frac{C^{t}_{i-1,j} - 2C^{t}_{i,j} + C^{t}_{i+1,j}}{\Delta x^2} \\ +\frac{\partial^2 C^t_{i,j}}{\partial y^2} &= \frac{C^{t}_{i,j-1} - 2C^{t}_{i,j} + C^{t}_{i,j+1}}{\Delta y^2} +\end{align} + +The ADI scheme is formally defined by the equations: + +#+NAME: eqn:genADI +\begin{equation} +\systeme{ + \displaystyle \frac{C^{t+1/2}_{i,j}-C^t_{i,j}}{\Delta t/2} = \displaystyle \alpha_x \frac{\partial^2 C^{t+1/2}_{i,j}}{\partial x^2} + \alpha_y \frac{\partial^2 C^{t}_{i,j}}{\partial y^2}, + \displaystyle \frac{C^{t+1}_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = \displaystyle \alpha_x \frac{\partial^2 C^{t+1/2}_{i,j}}{\partial x^2} + \alpha_y \frac{\partial^2 C^{t+1}_{i,j}}{\partial y^2} +} +\end{equation} + +\noindent The first of equations [[eqn:genADI]], which writes implicitly +the spatial derivatives in $x$ direction, after bringing the $\Delta t +/ 2$ terms on the right hand side and substituting $s_x=\frac{\alpha_x +\cdot \Delta t}{2 \Delta x^2}$ and $s_y=\frac{\alpha_y\cdot\Delta +t}{2\Delta y^2}$ reads: + +\begin{equation}\displaystyle +C^{t+1/2}_{i,j}-C^t_{i,j} = s_x (C^{t+1/2}_{i-1,j} - 2C^{t+1/2}_{i,j} + C^{t+1/2}_{i+1,j}) + s_y (C^{t}_{i,j-1} - 2C^{t}_{i,j} + C^{t}_{i,j+1}) +\end{equation} + +\noindent Separating the known terms (at time level $t$) on the left +hand side and the implicit terms (at time level $t+1/2$) on the right +hand side, we get: + +#+NAME: eqn:sweepX +\begin{equation}\displaystyle +-C^t_{i,j} - s_y (C^{t}_{i,j-1} - 2C^{t}_{i,j} + C^{t}_{i,j+1}) = - C^{t+1/2}_{i,j} + s_x (C^{t+1/2}_{i-1,j} - 2C^{t+1/2}_{i,j} + C^{t+1/2}_{i+1,j}) +\end{equation} + +\noindent Equation [[eqn:sweepX]] can be solved with a BTCS scheme since +it corresponds to a tridiagonal system of equations, and resolves the +$C^{t+1/2}$ at each inner grid cell. + +The second of equations [[eqn:genADI]] can be treated the same way and +yields: + +#+NAME: eqn:sweepY +\begin{equation}\displaystyle +-C^{t + 1/2}_{i,j} - s_x (C^{t + 1/2}_{i-1,j} - 2C^{t + 1/2}_{i,j} + C^{t + 1/2}_{i+1,j}) = - C^{t+1}_{i,j} + s_y (C^{t+1}_{i,j-1} - 2C^{t+1}_{i,j} + C^{t+1}_{i,j+1}) +\end{equation} + +This scheme only applies to inner cells, or else $\forall i,j \in [1, +n-1] \times [1, n-1]$. Following an analogous treatment as for the 1D +case, and noting $l_x$ and $l_y$ the constant left boundary values and +$r_x$ and $r_y$ the right ones for each direction $x$ and $y$, we can +modify equations [[eqn:sweepX]] for $i=0, j \in [1, n-1]$ + +#+NAME: eqn:boundXleft +\begin{equation}\displaystyle +-C^t_{0,j} - s_y (C^{t}_{0,j-1} - 2C^{t}_{0,j} + C^{t}_{0,j+1}) = - C^{t+1/2}_{0,j} + s_x (C^{t+1/2}_{1,j} - 3C^{t+1/2}_{0,j} + 2 l_x) +\end{equation} + +\noindent Similarly for $i=n, j \in [1, n-1]$: +#+NAME: eqn:boundXright +\begin{equation}\displaystyle +-C^t_{n,j} - s_y (C^{t}_{n,j-1} - 2C^{t}_{n,j} + C^{t}_{n,j+1}) = - C^{t+1/2}_{n,j} + s_x (C^{t+1/2}_{n-1,j} - 3C^{t+1/2}_{n,j} + 2 r_x) +\end{equation} + +\noindent For $i=j=0$: +#+NAME: eqn:bound00 +\begin{equation}\displaystyle +-C^t_{0,0} - s_y (C^{t}_{0,1} - 3C^{t}_{0,0} + 2l_y) = - C^{t+1/2}_{0,0} + s_x (C^{t+1/2}_{1,0} - 3C^{t+1/2}_{0,0} + 2 l_x) +\end{equation} + +Analogous expressions are readily derived for all possible +combinations of $i,j \in 0\times n$. In practice, wherever an index +$i$ or $j$ is $0$ or $n$, the centered spatial derivatives in $x$ or +$y$ directions must be substituted in relevant parts of the sweeping +equations \textbf{in both the implicit or the explicit sides} of +equations [[eqn:sweepX]] and [[eqn:sweepY]] by a term + +#+NAME: eqn:bound00 +\begin{equation}\displaystyle + s(C_{forw} - 3C + 2 bc) +\end{equation} +\noindent where $bc$ is the boundary condition in the given direction, +$s$ is either $s_x$ or $s_y$, and $C_{forw}$ indicates the contiguous +cell opposite to the boundary. Alternatively, noting the second +derivative operator as $\partial_{dir}^2$, we can write in compact +form: + +\begin{equation} +\systeme{ + \displaystyle \partial_x^2 C_{0,j} = 2l_x - 3C_{0,j} + C_{1,j} , + \displaystyle \partial_x^2 C_{n,j} = 2r_x - 3C_{n,j} + C_{n-1,j} , + \displaystyle \partial_y^2 C_{i,0} = 2l_y - 3C_{i,0} + C_{i,1} , + \displaystyle \partial_y^2 C_{i,n} = 2r_y - 3C_{i,n} + C_{i,n-1} +} +\end{equation} + + + +#+LATEX: \clearpage + +* Old stuff + +** Input - =c= $\rightarrow c$ - containing current concentrations at each grid cell for species @@ -18,7 +323,7 @@ - size: $N \times M$ - row-major -* Internals +** Internals - =A_matrix= $\rightarrow A$ - coefficient matrix for linear equation system implemented as sparse matrix @@ -34,7 +339,7 @@ - size: $(N+2) \cdot M$ - column-major (not relevant) -* Calculation for $\frac{1}{2}$ timestep +** Calculation for $\frac{1}{2}$ timestep ** Symbolic addressing of grid cells [[./grid.png]] @@ -77,9 +382,6 @@ s_x(i,j) & \text{else} - Adressing would look like this: $(i,j) = b(i \cdot (N+2) + j)$ - $\rightarrow$ for simplicity we will write $b(i,j)$ - - - *** Rules **** Ghost nodes @@ -96,11 +398,11 @@ c(i,N-1) & \text{else} *** Inlet -$p(i,j) = \frac{\Delta t}{2}\alpha(i,j)\frac{c(i-1,j) - 2\cdot c(i,j) + c(i+1,j)}{\Delta x^2}$[fn:1] +$p(i,j) = \frac{\Delta t}{2}\alpha(i,j)\frac{c(i-1,j) - 2\cdot c(i,j) + c(i+1,j)}{\Delta x^2}$ + +\noindent $p$ is called =t0_c= inside code $b(i,j) = \begin{cases} bc(i,j).\text{value} & \text{if } bc(i,N-1) = \text{constant} \\ -c(i,j)-p(i,j) & \text{else} \end{cases}$ - -[fn:1] $p$ is called =t0_c= inside code diff --git a/doc/grid_pqc.pdf b/doc/grid_pqc.pdf new file mode 100644 index 0000000..b970d18 Binary files /dev/null and b/doc/grid_pqc.pdf differ diff --git a/include/diffusion/BTCSDiffusion.hpp b/include/diffusion/BTCSDiffusion.hpp index ef4c4d6..e9eb868 100644 --- a/include/diffusion/BTCSDiffusion.hpp +++ b/include/diffusion/BTCSDiffusion.hpp @@ -135,16 +135,17 @@ private: const BCMatrixRowMajor &bc, double time_step, double dx) -> DMatrixRowMajor; - void fillMatrixFromRow(Eigen::SparseMatrix &A_matrix, - const DVectorRowMajor &alpha, - const BCVectorRowMajor &bc, int size, double dx, - double time_step); + static void fillMatrixFromRow(Eigen::SparseMatrix &A_matrix, + const DVectorRowMajor &alpha, + const BCVectorRowMajor &bc, int size, double dx, + double time_step); - void 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); + static void 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); void simulate3D(std::vector &c); inline static auto getBCFromFlux(Diffusion::boundary_condition bc, diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 8fdd020..240d76a 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -28,8 +28,6 @@ constexpr int BTCS_MAX_DEP_PER_CELL = 3; constexpr int BTCS_2D_DT_SIZE = 2; -constexpr double center_eq(double sx) { return -1. - 2. * sx; } - Diffusion::BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) { grid_cells.resize(dim, 1); @@ -104,7 +102,7 @@ void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c, b_vector.resize(size + 2); x_vector.resize(size + 2); - fillMatrixFromRow(A_matrix, alpha.row(0), bc.row(0), size, dx, time_step); + fillMatrixFromRow(A_matrix, alpha, bc, size, dx, time_step); fillVectorFromRow(b_vector, c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size, dx, time_step); @@ -143,17 +141,16 @@ void Diffusion::BTCSDiffusion::simulate2D( int n_rows = this->grid_cells[1]; int n_cols = this->grid_cells[0]; double dx = this->deltas[0]; - DMatrixRowMajor t0_c; double local_dt = this->time_step / BTCS_2D_DT_SIZE; - t0_c = calc_t0_c(c, alpha, bc, local_dt, dx); + DMatrixRowMajor t0_c = calc_t0_c(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), alpha.row(i), dx, local_dt, n_cols, - t0_c.row(i)); + simulate_base(input_field, bc.row(i + 1), alpha.row(i), dx, local_dt, + n_cols, t0_c.row(i)); c.row(i) << input_field; } @@ -165,8 +162,8 @@ void Diffusion::BTCSDiffusion::simulate2D( #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), alpha.col(i), dx, local_dt, n_rows, - t0_c.row(i)); + simulate_base(input_field, bc.col(i + 1), alpha.col(i), dx, local_dt, + n_rows, t0_c.row(i)); c.col(i) << input_field.transpose(); } } @@ -182,16 +179,22 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c, DMatrixRowMajor t0_c(n_rows, n_cols); - std::array y_values; + std::array y_values{}; // first, iterate over first row for (int j = 0; j < n_cols; j++) { - y_values[0] = getBCFromFlux(bc(0, j), c(0, j), alpha(0, j)); + boundary_condition tmp_bc = bc(0, j + 1); + + if (tmp_bc.type == Diffusion::BC_CLOSED){ + continue; + } + + 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) * - (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); + (2 * y_values[0] - 3 * y_values[1] + y_values[2]) / (dx * dx); } // then iterate over inlet @@ -212,12 +215,19 @@ 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; + } + y_values[0] = c(end - 1, j); y_values[1] = c(end, j); - y_values[2] = getBCFromFlux(bc(end, j), c(end, j), alpha(end, j)); + y_values[2] = getBCFromFlux(tmp_bc, c(end, j), alpha(end, j)); t0_c(end, j) = time_step * alpha(end, j) * - (y_values[0] - 2 * y_values[1] + y_values[2]) / (dx * dx); + (y_values[0] - 3 * y_values[1] + 2 * y_values[2]) / + (dx * dx); } return t0_c; @@ -227,8 +237,8 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( Eigen::SparseMatrix &A_matrix, const DVectorRowMajor &alpha, const BCVectorRowMajor &bc, int size, double dx, double time_step) { - Diffusion::boundary_condition left = bc[0]; - Diffusion::boundary_condition right = bc[size - 1]; + Diffusion::boundary_condition left = bc[1]; + Diffusion::boundary_condition right = bc[size]; bool left_constant = (left.type == Diffusion::BC_CONSTANT); bool right_constant = (right.type == Diffusion::BC_CONSTANT); @@ -239,27 +249,36 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( if (left_constant) { A_matrix.insert(1, 1) = 1; + } else { + double sx = (alpha[0] * time_step) / (dx * dx); + A_matrix.insert(1, 1) = -1. - 3. * sx; + A_matrix.insert(1, 0) = 2. * sx; + A_matrix.insert(1, 2) = sx; } - A_matrix.insert(A_size - 1, A_size - 1) = 1; - - if (right_constant) { - A_matrix.insert(A_size - 2, A_size - 2) = 1; - } - - for (int j = 1 + (int)left_constant, k = j - 1; - k < size - (int)right_constant; j++, k++) { + for (int j = 2, k = j - 1; k < size - 1; j++, k++) { double sx = (alpha[k] * time_step) / (dx * dx); - if (bc[k].type == Diffusion::BC_CONSTANT) { + if (bc[k + 1].type == Diffusion::BC_CONSTANT) { A_matrix.insert(j, j) = 1; continue; } - A_matrix.insert(j, j) = center_eq(sx); + A_matrix.insert(j, j) = -1. - 2. * sx; A_matrix.insert(j, (j - 1)) = sx; A_matrix.insert(j, (j + 1)) = sx; } + + if (right_constant) { + A_matrix.insert(A_size - 2, A_size - 2) = 1; + } else { + double sx = (alpha[size - 1] * time_step) / (dx * dx); + A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx; + A_matrix.insert(A_size - 2, A_size - 3) = sx; + A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx; + } + + A_matrix.insert(A_size - 1, A_size - 1) = 1; } void Diffusion::BTCSDiffusion::fillVectorFromRow( @@ -268,7 +287,7 @@ void Diffusion::BTCSDiffusion::fillVectorFromRow( const DVectorRowMajor &t0_c, int size, double dx, double time_step) { Diffusion::boundary_condition left = bc[0]; - Diffusion::boundary_condition right = bc[size - 1]; + Diffusion::boundary_condition right = bc[size + 1]; bool left_constant = (left.type == Diffusion::BC_CONSTANT); bool right_constant = (right.type == Diffusion::BC_CONSTANT); @@ -276,7 +295,7 @@ void Diffusion::BTCSDiffusion::fillVectorFromRow( int b_size = b_vector.size(); for (int j = 0; j < size; j++) { - boundary_condition tmp_bc = bc[j]; + boundary_condition tmp_bc = bc[j + 1]; if (tmp_bc.type == Diffusion::BC_CONSTANT) { b_vector[j + 1] = tmp_bc.value; @@ -287,15 +306,14 @@ void Diffusion::BTCSDiffusion::fillVectorFromRow( b_vector[j + 1] = -c[j] - t0_c_j; } - if (!left_constant) { - // this is not correct currently.We will fix this when we are able to define - // FLUX boundary conditions - b_vector[0] = getBCFromFlux(left, c[0], alpha[0]); - } + // this is not correct currently.We will fix this when we are able to define + // FLUX boundary conditions + b_vector[0] = + (left_constant ? left.value : getBCFromFlux(left, c[0], alpha[0])); - if (!right_constant) { - b_vector[b_size - 1] = getBCFromFlux(right, c[size - 1], alpha[size - 1]); - } + b_vector[b_size - 1] = + (right_constant ? right.value + : getBCFromFlux(right, c[size - 1], alpha[size - 1])); } void Diffusion::BTCSDiffusion::setTimestep(double time_step) { @@ -312,7 +330,7 @@ auto Diffusion::BTCSDiffusion::simulate(double *c, double *alpha, if (this->grid_dim == 1) { Eigen::Map c_in(c, this->grid_cells[0]); Eigen::Map alpha_in(alpha, this->grid_cells[0]); - Eigen::Map bc_in(bc, this->grid_cells[0]); + Eigen::Map bc_in(bc, this->grid_cells[0] + 2); simulate1D(c_in, alpha_in, bc_in); } @@ -323,8 +341,8 @@ auto Diffusion::BTCSDiffusion::simulate(double *c, double *alpha, Eigen::Map alpha_in(alpha, this->grid_cells[1], this->grid_cells[0]); - Eigen::Map bc_in(bc, this->grid_cells[1], - this->grid_cells[0]); + Eigen::Map bc_in(bc, this->grid_cells[1] + 2, + this->grid_cells[0] + 2); simulate2D(c_in, alpha_in, bc_in); }