From 9faa4e79bbe8f62a731076dc831dfc96c1102f08 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 6 Apr 2022 09:50:08 +0200 Subject: [PATCH 01/36] Prepare setup of matrix A for new equation --- src/BTCSDiffusion.cpp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 8fdd020..7277874 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); @@ -239,16 +237,25 @@ 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. - 2. * sx; + A_matrix.insert(1, 0) = 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; + } else { + double sx = (alpha[size - 1] * time_step) / (dx * dx); + A_matrix.insert(A_size - 2, A_size - 2) = -1. - 2. * sx; + A_matrix.insert(A_size - 2, A_size - 3) = sx; + A_matrix.insert(A_size - 2, A_size - 1) = sx; } - 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) { @@ -256,7 +263,7 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( 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; } From ad6e1ad616a59d23c49f08504aae66035b142e81 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 19 Apr 2022 10:39:37 +0200 Subject: [PATCH 02/36] Allow boundary conditions in ghost nodes too --- src/BTCSDiffusion.cpp | 42 +++++++++++++++++++++--------------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 7277874..1271156 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -102,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); @@ -225,8 +225,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); @@ -244,7 +244,18 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( A_matrix.insert(1, 2) = sx; } - A_matrix.insert(A_size - 1, A_size - 1) = 1; + for (int j = 2, k = j - 1; k < size - 1; j++, k++) { + double sx = (alpha[k] * time_step) / (dx * dx); + + if (bc[k + 1].type == Diffusion::BC_CONSTANT) { + A_matrix.insert(j, j) = 1; + continue; + } + + 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; @@ -255,18 +266,7 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( A_matrix.insert(A_size - 2, A_size - 1) = sx; } - 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) { - A_matrix.insert(j, j) = 1; - continue; - } - - A_matrix.insert(j, j) = -1. - 2. * sx; - A_matrix.insert(j, (j - 1)) = sx; - A_matrix.insert(j, (j + 1)) = sx; - } + A_matrix.insert(A_size - 1, A_size - 1) = 1; } void Diffusion::BTCSDiffusion::fillVectorFromRow( @@ -275,7 +275,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); @@ -283,7 +283,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; @@ -319,7 +319,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); } @@ -330,8 +330,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); } From 7b36225bd651b48a9218cb624ec056924a9170c9 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 19 Apr 2022 10:49:18 +0200 Subject: [PATCH 03/36] Update app to new API --- app/main_1D.cpp | 6 +++--- app/main_2D.cpp | 10 +++++----- app/main_2D_mdl.cpp | 22 +++++++++++----------- 3 files changed, 19 insertions(+), 19 deletions(-) diff --git a/app/main_1D.cpp b/app/main_1D.cpp index b2423f2..e908039 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[1] = {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_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; } From 777d75baa557ee554f59d1284eb2c91f9003307c Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 09:34:32 +0200 Subject: [PATCH 04/36] Checkout files from age --- doc/ADI_scheme.org | 94 +++++++++++++++++++++++++++++++++++++++++++++ doc/grid_pqc.pdf | Bin 0 -> 5809 bytes 2 files changed, 94 insertions(+) create mode 100644 doc/grid_pqc.pdf diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 6a86caf..8a3fd92 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -1,6 +1,7 @@ #+TITLE: Adi 2D Scheme #+LaTeX_CLASS_OPTIONS: [a4paper,10pt] #+LATEX_HEADER: \usepackage{fullpage} +#+LATEX_HEADER: \usepackage{amsmath} #+OPTIONS: toc:nil * Input @@ -104,3 +105,96 @@ bc(i,j).\text{value} & \text{if } bc(i,N-1) = \text{constant} \\ \end{cases}$ [fn:1] $p$ is called =t0_c= inside code + +** Finite differences with nodes as cells' centres + +*** The explicit FTCS scheme as in PHREEQC + +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 discretize it following a Forward Time, Centered Space finite +difference scheme where the nodes correspond to the centers of a grid +such as: + +[[./grid_pqc.pdf]] + +The left boundary is defined on $x=0$ while the center of the first +cell is in $x=dx/2$, with $dx=L/n$. + +We discretize [[eqn:1]] as following, for each index i in 1, \dots, n-1 +and assuming constant $\alpha$: + +#+NAME: eqn:2 +\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} + +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 +boundaries. 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^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}- +\frac{C^{j+1}_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +#+NAME: eqn:4 +\begin{align}\displaystyle +C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}-C^{j+1}_{0}- 2 C^{j+1}_{0}+2l \right) \nonumber \\ + & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) +\end{align} + +Which leads to the following term applicable to our approach: + +#+NAME: eqn:5 +\begin{equation}\displaystyle + -C^j_0} = \left[\frac{\alpha \cdot \Delta t}{\Delta x^2}\right] \cdot \left(C^{j+1}_1 + 2l \right) + \left[ -1 - 3 \cdot \frac{\alpha \cdot \Delta t}{\Delta x^2} \right] \cdot C^{j+1}_0 +\end{equation} + +In case of constant right boundary, the finite difference of point +$C_n$ - calling $r$ the right boundary value - is: + +#+NAME: eqn:6 +\begin{equation}\displaystyle +\frac{C_n^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}- +\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +Which, developed, gives +#+NAME: eqn:7 +\begin{align}\displaystyle +C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\ + & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{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:8 +\begin{equation}\displaystyle +C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n) +\end{equation} + + + +A similar treatment can be applied to the BTCS implicit scheme. diff --git a/doc/grid_pqc.pdf b/doc/grid_pqc.pdf new file mode 100644 index 0000000000000000000000000000000000000000..b970d186492f1e3f6c90268bc2a56447fc4cabe8 GIT binary patch literal 5809 zcma)A2{@Ep*p?d0WG_2;?IFw=mSt47D>`}JJ zk|k@1lBHyee}?M&zU%Az|9{?@_q^vh&pGEk&vQR>UH2_$bVgGgDj^FNto&TE4VDHV z0M6ke7=ZvJ^-)+Cf-3+gnoPg|07z=0-3Ta8V(W$^pw6HiaZV_(k`fqC@I)aoV4t*W zrUr4xpfu0tEH;MZV=gI-Q%Y!4p+S`P*mGgGKUW0~A4Umot$9eh{kFFI76#mySQNMj<0=Z7Si1tmoY-}345?G@#Fm~24nYRT zF(A&1F>r@Z)N5LIGQ8Cfn~O$f6gI?O`b1vwKuAzdO#7L1paWcWv6WBEeq9>4_d!Z< z#_YC9l%i09Cg{;pN_0`PZauSml1j{lyFACjq{{7JVM4U{^(j}j@}`;D-rGlPJ;l`T z&6YtrVv`#3G@o{sJ{}o0D{W1!PU|OdA!wa1mb06_6}s=oCUWG)xL_KRLf|;)Z5H~0 z@%M@^_Xs!Qa$aY$GQNH|yyk$aD7oqR@<}tNpgjMFwjA#V~o=}`6Zo?mZ(qg?9P?1-tb$J)11onb_B_e!o@<_JA~c#i@duPEom zv7AlH&Y+RIC3DZRz~!M#j$J1%1PM(o#;N&j?V8+9fzR%h`R8X|$gW~9j`{);xeexD zs=9w|kTN|F>udLhP9w+R@B%kqIsdx*VzyE5W$Oh!&n|Hh#f#uT@`==l8#8?oro9W zz`u}MOfZ5;rCsO5p&6G3* zp!)-wXiq!=kpB7nzoS>N^EtQ?avU`A{HzDdoKWTIFgi?@9F@sfGjDS=K}XlWx$e1+ zKRqgF@kSdH~Bu*ivb|L;gAqymwVnE-<7onMW&{pO4kj>xT{<->yMUt zcDEDZ$S?+xoSl)j^Ya}2x}OIg=TQ^z2>XE~g%hg&k4jEAr?xaZ z$Eh^_<4$c~U=3ZPeup5Y8)N``olEX*L&nU=pn6D^O@`5m!r9)Ef`Y<{$L!2IRfdjN z;;^K*uWT7B+gR}J9aZ>`Q{w?|cDJq+KDbvmS34&^H~K|kZo8}Y`!zejT${{(D7YuG zzQ!`HEt0I$H7+y0W;u|J+&q+lA*0|~)N&kosk0BA?-em*RN3@{nT%n<*ida}nol?_ z`URXLMYMRzNtyC#bssIxDW{M#Z_+R9vR%Z~&)pNv1-v)yT$N;_dT_f_xmo-*ZO_D6={ zz#~2;KEYErZ)!!RSsvq*3LHuZa$-nPpH1*>zR@9kv`wu%vgR25lnzUxpp(c7-|I+O zE5oaik}LdTk;>G9NGEY>Va7A+d^)nV`X%v2GyJR2kGvW2mu5^?**DX!s+J^5w$pkh zj`CHvTg=dHrj80epAqJZSMto*iP)j}QQbrq%&e^fKGinC7#{*3swXQ83xZHlG%>J& zn^X>YLAX9O9J-YmtkuBywm#@sOgI$6l7{4piI_u@>xZ)1vvh>ZLs(~nI6ERhWehAq z3O<3tr@{rxjz6lyEC-?{rbITx`2+2u`(NJr5&bnZQW%y_qZsNq%qvH+Vmk;ky9ru|J%TNha zs-fQt%3yCKYl5^{LYWrFImbq}AF8>I0g4$`w-e_07OkcGZsJhi-wX;sO~x6;?t?UL zyHyd+U-0ucyc-i|zt)wU?0pwC5P*I7Mo?ohL2zXG^+rCQupqyzUZqemf&S%-Zs*xe zp|FKnfztwy^0I?R9s0W*t^vYw6`iWaxkYB=3Yzg${GmuzAtS`-g-UNq**LgM-vc_` z(7E@fY%~Ra>8F*T=Sz(Z&U=Q@Wt(qH61a!wp5EGyzB>B#K|x+hIVNCbm)@B{T5J9~ z2giW1(cZ{Iz6H*g3pHD<3k%`lm-{)o*1xcdCilhSc0ZQ9NR5xT{#I6vwauATrR@6N z15$Jyv+Iz*(7XJ=gT>3UuH#n1>gZ9i&_I&!@&Uaef`kZ)od zpG!-?DQbx(+LgPBu)}PiHI8HL(vT}Du#He7hf6)OYXia)1o|6$G#}^R7esYDQ^A;} zo{b6-vgPw<4%@8ND70-sAF{C%->{W5oxzg6|@G-wR?DGFsYL< z^QHn@2$3fpNzrQ#@G0_c>p+=z#>&ysh`QBIcSAPKZ3!)c!Ay>sO!d9~lvtgc#MxO+YY_g5DL407*_)gys-lGGcdf6Ej|(m3E!P0d zLRv=_Xp!&e+P_jCEqcRz3wP>N>Ch|Nyu+nO4Vmzo7xGk7?AySwLQI^|UHXDm^@+Pa z3K|(`I$Q!+wF38RP4Q@Sw9finSo(07Lg-`^ z8%OVmNNBT-qxX!Nz+E+`@4{vBmo7gkPR5`Q*bQYa;}c;Q_1rO!-f zu~lZ;HOi)*x`_T!3I7pL=3pr=j}AAp2&iG{bOcAKaA_@@R`Bllyz^fN|_AfYo0@ny~kfXN^He9>NnUXyb`mr#(-Wt4lT)#uX?Rq zuiX zq{v38*`9w#i(N~z(>>Qh_>~Vamd{Y{#3p$r_v_m5=y0)1BOHL8%ICti;8ny0C4R)nq&&QgaT!B?`43v}h zT*2xdY@$7M5%(2O0}+NBY_at>H{_Ip&Q6(bNGpYOPvPajTGMP?(Oq(MFQ?b!4hc^k zm7}(6Zr@<4i!pTp&EO0SkVJ2(U3bI_PU2N`9A$tHcyF4fX($o zY>`C}ayZ*{W40Smw&+C&-Jl}^TYp5qSwLS?JkMTq!6*)WCzhh4UcaN!P1=EOHk*14 z!qaIP7*@j=UB;$YcKmf2msn45QVq52G|R)rV#^yN?g|6;;sZ!U8GC*iB&=$J&aFje z^g3?x#Ax)kyYh+p>J^|aBx2bX5u>$y(rfC94^7-O^;I6_dioWCdyRcQ6x~zuh&a|| zKCel^Ppz+gK&4Y!8_@L$VINxMmhNTz4teQxvd>}Ml=!n7c9YyY9ABc&>}5MV>DKK; zB)xnU&d(u{HP+*XkhuRc=OfQNr$Z_V!{~4Q-0XcFb(U%x&i>Y-Q*Ne`?9K9O2``g! zYjHxKanw>Tw?bS&EG90^sLRT(ut-ZOrrLs_?>KmAO|}*9p^=a<-xloJf5*c1<&%Dj zbeidG$HA8m7RoEQ2414(nUXwmpNy?sWG$U9v;WYmw_7dgHMpddkHrnY9W6mDZEqDh zS7x5Pb*epDyu{4B^6>z3LB&WF+nSkzx79*w&22xu#X~6dSs6^ldB|50=5X#Yc3e@N zs}8kCM&0@NnD4GWcLqlcK}^0G&dT(OLsBoy^Tzm2Jw3LF#}nwc z&cTCeTIAW*!kc(CsY_Sk2FhYj?Un%L;L~Ex8==j9 z`K?cEvp#u$xmP=EUP3QA37!^6&F;M3UpQQBUB!Lxk$3MZYK1o_tvX=lVRBON)GAvP z?sK8rMbkmdh{^cSrV4_2$T&pfw1Q+o$?CV6wAm zNya?u+x6;+te|Fzucm$cZqB=@Z|{6<(V;V@-xy9yw$AsvfAX5(V70kfK+4Ax6joBS z)n!YdkGYyqM95ghQ;$B0{@v7g_O!Az_>;Vk2q}#WxZ7SV=S(7ph1U{p+`mjC-9+@= z^u{F#!PXu9%I%y;BZF*xm$tSxBj>v@o0{JgUYxDP>p5PzTzjFyirdyIZ#-wX&UdX! zW^XtQd%;*NyLO)WBd@&VYoWa`rqx?VBXRN-aW!RxXIYn@*w z>|rYd9JqJJL^rnHzo6ZdtdRaj|8Ewt&kcUk2pB{P`io*grT)MpNF?w7$Np^-y0IY; zn$|DAV4kPypiuJS(D=NiMr43=6=h)aV{QI3yo$@~_^HA~FS5Ohhds zLx2qM_2w9Jk|^6&@=ja84(v|Ew5p4yQb@-Nd>lJ zk^V@q($xl2+VgETINPv~V!^}h2*(fHSzRo<=H<`n`Fc9JXx$(=m#*v4P$|$Wv96+e z5T3EdN#c5ux)V72Hrfe%OWAW*J8iM;cx0_m9Q;#O)Ct3fzF-@uMUd44n$wi)cq;#> zplftNRxcvA*cV;0lAkBN;oM86*-KjU770}y6K2w6>00+-S(xSGR)0q3v5Q19k42drr znCqMBiy5LZ4wvyJI1CaiZh~^T?1ucyQ4@(lyZHhq4}5@qCy|3?+>kE#UqqU?q?#Je z2e1ajp)zuSxRe|efI=V;;tdW1Z1*SV5Rh(YM^&tg8w!BzTU7CmC?egLk%R4bN$*MO zi^HI@`v!HSyEY2#;!47+zri-HlI_YwM+Qb5uVNa!XB z4}XQ_V7-4(-hGoA5|7%KA^y+!kfa24SEMH~2m3Sikw5!Tse>2%XY>ImDL)6Lqk$!A zWoWDmAZdZds$%izUxxjKEYMB_S3EJAQgG5){`-=Wl>^|!F#r1f$rCO^+V;OMnf>bj z?{OH6gtmh_X&=%>YNe&6f48L3U++MACz8*v_F%@r{J(4@?|oblQ-V8(MU$jMVE_Nd z({Q>?M&_|Rd$yKkz{m`&$|DZC~vTWd}WvOkjY}6fV@kZ z`TL<<(+m7dqBnlC0<;1t6ARF3f2Xw4rz<_PEeL;e5KzKHb? zm4W@P|6gj*UqH7sK{+2_^FX?lmD}(ADJJ_05ZaHJ1^ELI0@hGK3XnReUK1RSScboq zI7uC>GY)_rD8YY^rlY<6#$Lwtpv0h{>p|&3K{wkuv?0qOmW3*6^GuPI Date: Wed, 20 Apr 2022 09:55:38 +0200 Subject: [PATCH 05/36] Apply new scheme to model (only 1D) --- src/BTCSDiffusion.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 1271156..4df37a9 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -239,8 +239,8 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( A_matrix.insert(1, 1) = 1; } else { double sx = (alpha[0] * time_step) / (dx * dx); - A_matrix.insert(1, 1) = -1. - 2. * sx; - A_matrix.insert(1, 0) = sx; + A_matrix.insert(1, 1) = -1. - 3. * sx; + A_matrix.insert(1, 0) = 2.*sx; A_matrix.insert(1, 2) = sx; } @@ -261,9 +261,9 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( 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. - 2. * sx; + 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) = sx; + A_matrix.insert(A_size - 2, A_size - 1) = 2.*sx; } A_matrix.insert(A_size - 1, A_size - 1) = 1; From 899e00bec1e5d0ee1139b6d814758efab2d691d0 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 10:06:38 +0200 Subject: [PATCH 06/36] Update applicable equation --- doc/ADI_scheme.org | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 8a3fd92..557fe34 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -163,11 +163,17 @@ C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) \end{align} -Which leads to the following term applicable to our approach: +Now we define variable $s_x$ as following: + +\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: #+NAME: eqn:5 \begin{equation}\displaystyle - -C^j_0} = \left[\frac{\alpha \cdot \Delta t}{\Delta x^2}\right] \cdot \left(C^{j+1}_1 + 2l \right) + \left[ -1 - 3 \cdot \frac{\alpha \cdot \Delta t}{\Delta x^2} \right] \cdot C^{j+1}_0 + -C^j_0} = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} In case of constant right boundary, the finite difference of point From 0ef5568e102a6377e2e82ea10fba50b1c92cf67a Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 10:09:43 +0200 Subject: [PATCH 07/36] Fix broken reference --- doc/ADI_scheme.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 557fe34..3cabbf1 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -192,7 +192,7 @@ C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{n-1} \right) \end{align} -If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:5]] +If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:6]] becomes zero and we are left with: From b251204bf5a3d915d3918d3ad97b34d201fac905 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 11:27:30 +0200 Subject: [PATCH 08/36] Apply bc to ghost node --- app/main_1D.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/app/main_1D.cpp b/app/main_1D.cpp index e908039..2fb7c21 100644 --- a/app/main_1D.cpp +++ b/app/main_1D.cpp @@ -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[1] = {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)); From f1ca3ff2eaccb02c2f4100db22e37976c97fd788 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 11:30:59 +0200 Subject: [PATCH 09/36] Revert to 'age' state and append with implicit BTCS --- doc/ADI_scheme.org | 96 ++++++++++++++++++++++++++++++---------------- 1 file changed, 64 insertions(+), 32 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 3cabbf1..09158ed 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -132,7 +132,7 @@ and assuming constant $\alpha$: #+NAME: eqn:2 \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} + \frac{C_i^{j+1} -C_i^{j}}{\Delta t} = \alpha\frac{\frac{C^j_{i+1}-C^j_{i}}{\Delta x}-\frac{C^j_{i}-C^j_{i-1}}{\Delta x}}{\Delta x} \end{equation} In practice, we evaluate the first derivatives of $C$ w.r.t. $x$ on @@ -150,6 +150,69 @@ 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^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^j_{1}-C^j_{0}}{\Delta x}- +\frac{C^j_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +#+NAME: eqn:4 +\begin{align}\displaystyle +C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^j_{1}-C^j_{0}- 2 C^j_{0}+2l \right) \nonumber \\ + & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^j_{1}- 3 C^j_{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^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}- +\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +Which, developed, gives +#+NAME: eqn:6 +\begin{align}\displaystyle +C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\ + & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{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^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n) +\end{equation} + + + +A similar treatment can be applied to the BTCS implicit scheme. + +*** implicit BTCS + +\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} + +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 +boundaries. 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: + \begin{equation}\displaystyle \frac{C_0^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}- \frac{C^{j+1}_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} @@ -157,7 +220,6 @@ calling $l$ the numerical value of a constant boundary condition: This expression, once developed, yields: -#+NAME: eqn:4 \begin{align}\displaystyle C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}-C^{j+1}_{0}- 2 C^{j+1}_{0}+2l \right) \nonumber \\ & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) @@ -171,36 +233,6 @@ Now we define variable $s_x$ as following: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: -#+NAME: eqn:5 \begin{equation}\displaystyle -C^j_0} = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} - -In case of constant right boundary, the finite difference of point -$C_n$ - calling $r$ the right boundary value - is: - -#+NAME: eqn:6 -\begin{equation}\displaystyle -\frac{C_n^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}- -\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x} -\end{equation} - -Which, developed, gives -#+NAME: eqn:7 -\begin{align}\displaystyle -C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\ - & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{n-1} \right) -\end{align} - -If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:6]] -becomes zero and we are left with: - - -#+NAME: eqn:8 -\begin{equation}\displaystyle -C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n) -\end{equation} - - - -A similar treatment can be applied to the BTCS implicit scheme. From ed19c5a604dc6b22a71dd8941c17767192f96aac Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Thu, 21 Apr 2022 12:51:34 +0200 Subject: [PATCH 10/36] Fix unnecessary bracket --- doc/ADI_scheme.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 09158ed..d47f3cb 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -234,5 +234,5 @@ Now we define variable $s_x$ as following: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: \begin{equation}\displaystyle - -C^j_0} = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + -C^j_0 = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} From 9da59af3f77a3b60dd38a5ba79bde518cf16e59e Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Fri, 22 Apr 2022 09:49:53 +0200 Subject: [PATCH 11/36] Update documentation of implicit BTCS --- doc/ADI_scheme.org | 41 +++++++++++++++++++++++++++-------------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index d47f3cb..7a0b0c7 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -195,23 +195,36 @@ A similar treatment can be applied to the BTCS implicit scheme. *** implicit BTCS -\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} +First, we define the Backward time difference: + +\begin{equation} + \frac{\partial C }{\partial t} = \frac{C^j_i - C^{j-1}_i}{\Delta t} \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$. +Second the spatial derivative approximation: -This discretization works for all internal cells, but not for the -boundaries. To properly treat them, we need to account for the -discrepancy in the discretization. +\begin{equation} + \frac{\partial^2 C }{\partial t} = \frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{i-1}}{\Delta x}}{\Delta x} +\end{equation} -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: +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} -C_i^{j-1}}{\Delta t} = \alpha\frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{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: + +\begin{align}\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} \nonumber \\ + & = \alpha\frac{C^{j+1}_{i-1} - 2C^{j+1}_{i} + C^{j+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^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}- @@ -225,7 +238,7 @@ C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) \end{align} -Now we define variable $s_x$ as following: +Now we define variable $s_x$ as followed: \begin{equation} s_x = \frac{\alpha \cdot \Delta t}{\Delta x^2} From 82da351bba609a70571700d466a4cf1174fcfa82 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Fri, 22 Apr 2022 13:55:50 +0200 Subject: [PATCH 12/36] MDL: some restructuring in ADI_scheme.org --- doc/ADI_scheme.org | 259 ++++++++++++++++++++++++--------------------- 1 file changed, 136 insertions(+), 123 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 7a0b0c7..a3497aa 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -1,114 +1,11 @@ -#+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} #+OPTIONS: toc:nil -* Input -- =c= $\rightarrow c$ - - containing current concentrations at each grid cell for species - - size: $N \times M$ - - row-major -- =alpha= $\rightarrow \alpha$ - - diffusion coefficient for both directions (x and y) - - size: $N \times M$ - - row-major -- =boundary_condition= $\rightarrow bc$ - - Defines closed or constant boundary condition for each grid cell - - size: $N \times M$ - - row-major - -* Internals - -- =A_matrix= $\rightarrow A$ - - coefficient matrix for linear equation system implemented as sparse matrix - - size: $((N+2)\cdot M) \times ((N+2)\cdot M)$ (including ghost zones in x direction) - - column-major (not relevant) - -- =b_vector= $\rightarrow b$ - - right hand side of the linear equation system - - size: $(N+2) \cdot M$ - - column-major (not relevant) -- =x_vector= $\rightarrow x$ - - solutions of the linear equation system - - size: $(N+2) \cdot M$ - - column-major (not relevant) - -* Calculation for $\frac{1}{2}$ timestep - -** Symbolic addressing of grid cells -[[./grid.png]] - -** Filling of matrix $A$ - -- row-wise iterating with $i$ over =c= and =\alpha= matrix respectively -- addressing each element of a row with $j$ -- matrix $A$ also containing $+2$ ghost nodes for each row of input matrix $\alpha$ - - $\rightarrow offset = N+2$ - - addressing each object $(i,j)$ in matrix $A$ with $(offset \cdot i + j, offset \cdot i + j)$ - -*** Rules - -$s_x(i,j) = \frac{\alpha(i,j)*\frac{t}{2}}{\Delta x^2}$ where $x$ defining the domain size in x direction. - -For the sake of simplicity we assume that each row of the $A$ matrix is addressed correctly with the given offset. - -**** Ghost nodes - -$A(i,-1) = 1$ - -$A(i,N) = 1$ - -**** Inlet - -$A(i,j) = \begin{cases} -1 & \text{if } bc(i,j) = \text{constant} \\ --1-2*s_x(i,j) & \text{else} -\end{cases}$ - -$A(i,j\pm 1) = \begin{cases} -0 & \text{if } bc(i,j) = \text{constant} \\ -s_x(i,j) & \text{else} -\end{cases}$ - -** Filling of vector $b$ - -- each elements assign a concrete value to the according value of the row of matrix $A$ -- 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 - -$b(i,-1) = \begin{cases} -0 & \text{if } bc(i,0) = \text{constant} \\ -c(i,0) & \text{else} -\end{cases}$ - -$b(i,N) = \begin{cases} -0 & \text{if } bc(i,N-1) = \text{constant} \\ -c(i,N-1) & \text{else} -\end{cases}$ - -*** 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] - -$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 - -** Finite differences with nodes as cells' centres - -*** The explicit FTCS scheme as in PHREEQC +* Finite differences with nodes as cells' centres The 1D diffusion equation is: @@ -118,17 +15,22 @@ The 1D diffusion equation is: & = \alpha \frac{\partial^2 C}{\partial x^2} \end{align} -We discretize it following a Forward Time, Centered Space finite -difference scheme where the nodes correspond to the centers of a grid -such as: +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 is in $x=dx/2$, with $dx=L/n$. +cell - which are the points constituting the finite difference nodes - +is in $x=dx/2$, with $dx=L/n$. -We discretize [[eqn:1]] as following, for each index i in 1, \dots, n-1 -and assuming constant $\alpha$: + +** 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 @@ -142,8 +44,8 @@ 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 -boundaries. To properly treat them, we need to account for the -discrepancy in the discretization. +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, @@ -193,29 +95,30 @@ C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} A similar treatment can be applied to the BTCS implicit scheme. -*** implicit BTCS +** Implicit BTCS scheme -First, we define the Backward time difference: +First, we define the Backward Time difference: \begin{equation} - \frac{\partial C }{\partial t} = \frac{C^j_i - C^{j-1}_i}{\Delta t} + \frac{\partial C^{j+1} }{\partial t} = \frac{C^{j+1}_i - C^{j}_i}{\Delta t} \end{equation} -Second the spatial derivative approximation: +Second the spatial derivative approximation, evaluated at time level $j+1$: \begin{equation} - \frac{\partial^2 C }{\partial t} = \frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{i-1}}{\Delta x}}{\Delta x} + \frac{\partial^2 C^{j+1} }{\partial x^2} = \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} 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} -C_i^{j-1}}{\Delta t} = \alpha\frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{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: +# \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: \begin{align}\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} \nonumber \\ @@ -249,3 +152,113 @@ Substituting with the new variable $s_x$ and reordering of terms leads to the eq \begin{equation}\displaystyle -C^j_0 = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} + +*TODO* +- Right boundary +- Tridiagonal matrix filling + + + + +#+LATEX: \clearpage + +* Old stuff + +** Input + +- =c= $\rightarrow c$ + - containing current concentrations at each grid cell for species + - size: $N \times M$ + - row-major +- =alpha= $\rightarrow \alpha$ + - diffusion coefficient for both directions (x and y) + - size: $N \times M$ + - row-major +- =boundary_condition= $\rightarrow bc$ + - Defines closed or constant boundary condition for each grid cell + - size: $N \times M$ + - row-major + +** Internals + +- =A_matrix= $\rightarrow A$ + - coefficient matrix for linear equation system implemented as sparse matrix + - size: $((N+2)\cdot M) \times ((N+2)\cdot M)$ (including ghost zones in x direction) + - column-major (not relevant) + +- =b_vector= $\rightarrow b$ + - right hand side of the linear equation system + - size: $(N+2) \cdot M$ + - column-major (not relevant) +- =x_vector= $\rightarrow x$ + - solutions of the linear equation system + - size: $(N+2) \cdot M$ + - column-major (not relevant) + +** Calculation for $\frac{1}{2}$ timestep + +** Symbolic addressing of grid cells +[[./grid.png]] + +** Filling of matrix $A$ + +- row-wise iterating with $i$ over =c= and =\alpha= matrix respectively +- addressing each element of a row with $j$ +- matrix $A$ also containing $+2$ ghost nodes for each row of input matrix $\alpha$ + - $\rightarrow offset = N+2$ + - addressing each object $(i,j)$ in matrix $A$ with $(offset \cdot i + j, offset \cdot i + j)$ + +*** Rules + +$s_x(i,j) = \frac{\alpha(i,j)*\frac{t}{2}}{\Delta x^2}$ where $x$ defining the domain size in x direction. + +For the sake of simplicity we assume that each row of the $A$ matrix is addressed correctly with the given offset. + +**** Ghost nodes + +$A(i,-1) = 1$ + +$A(i,N) = 1$ + +**** Inlet + +$A(i,j) = \begin{cases} +1 & \text{if } bc(i,j) = \text{constant} \\ +-1-2*s_x(i,j) & \text{else} +\end{cases}$ + +$A(i,j\pm 1) = \begin{cases} +0 & \text{if } bc(i,j) = \text{constant} \\ +s_x(i,j) & \text{else} +\end{cases}$ + +** Filling of vector $b$ + +- each elements assign a concrete value to the according value of the row of matrix $A$ +- 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 + +$b(i,-1) = \begin{cases} +0 & \text{if } bc(i,0) = \text{constant} \\ +c(i,0) & \text{else} +\end{cases}$ + +$b(i,N) = \begin{cases} +0 & \text{if } bc(i,N-1) = \text{constant} \\ +c(i,N-1) & \text{else} +\end{cases}$ + +*** 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}$ + +\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}$ From 39ce2a62c3a2aeea0ecf461e6ef702a46a59d361 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 25 Apr 2022 10:52:23 +0200 Subject: [PATCH 13/36] Reorder equation for better readability --- doc/ADI_scheme.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index a3497aa..d363194 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -150,7 +150,7 @@ Now we define variable $s_x$ as followed: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: \begin{equation}\displaystyle - -C^j_0 = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + -C^j_0 = (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + s_x \cdot C^{j+1}_1 \end{equation} *TODO* From 9a5b7c1101ede7240e0a303fd44433a021d180ad Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 25 Apr 2022 10:50:06 +0200 Subject: [PATCH 14/36] Added right boundary --- doc/ADI_scheme.org | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index d363194..391d3a6 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -153,8 +153,27 @@ Substituting with the new variable $s_x$ and reordering of terms leads to the eq -C^j_0 = (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + s_x \cdot C^{j+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^{j+1} -C_n^{j}}{\Delta t} = \alpha\frac{\frac{r-C^{j+1}_{n}}{\frac{\Delta x}{2}}- +\frac{C^{j+1}_{n}-C^{j+1}_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +\begin{align}\displaystyle +C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 2C^{j+1}_{n} - C^{j+1}_{n} + C^{j+1}_{n-1} \right) \nonumber \\ + & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 3C^{j+1}_{n} + C^{j+1}_{n-1} \right) +\end{align} + +Now rearrange terms and substituting with $s_x$ leads to: + +\begin{equation}\displaystyle + -C^j_n = s_x \cdot C^{j+1}_{n-1} + (-1 - 3s_x) \cdot C^{j+1}_n + (2s_x) \cdot r +\end{equation} + *TODO* -- Right boundary - Tridiagonal matrix filling From 6de575ed8f96fb32ac386aaafbc008fc751ef733 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 25 Apr 2022 12:31:27 +0200 Subject: [PATCH 15/36] Added 2D ADI scheme --- doc/ADI_scheme.org | 48 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 391d3a6..3f3e303 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -105,6 +105,7 @@ First, we define the Backward Time difference: Second the spatial derivative approximation, evaluated at time level $j+1$: +#+NAME: eqn:secondderivative \begin{equation} \frac{\partial^2 C^{j+1} }{\partial x^2} = \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} @@ -120,6 +121,7 @@ equations given above leads to the following 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^{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} \nonumber \\ & = \alpha\frac{C^{j+1}_{i-1} - 2C^{j+1}_{i} + C^{j+1}_{i+1}}{\Delta x^2} @@ -176,8 +178,54 @@ Now rearrange terms and substituting with $s_x$ leads to: *TODO* - Tridiagonal matrix filling +** 2D ADI using BTCS scheme +In the previous sections we described the usage of FTCS and BTCS on 1D grids. To +make usage of the BTCS scheme we are following the alternating-direction +implicit method. Therefore we make use of second difference operator defined in +equation [[eqn:secondderivative]] in both $x$ and $y$ direction for half the time +step $\Delta t$. We are denoting the numerator of equation [[eqn:1DBTCS]] as +$\delta^2_x C^n_{ij}$ and $\delta^2_y C^n_{ij}$ respectively with $i$ the +position in $x$, $j$ the position in $y$ direction and $n$ as the current time +step: +\begin{align}\displaystyle +\delta^2_x C^n_{ij} &= C^{n}_{i-1,j} - 2C^{n}_{i,j} + C^{n}_{i+1,j} \nonumber \\ +\delta^2_y C^n_{ij} &= C^{n}_{i,j-1} - 2C^{n}_{i,j} + C^{n}_{i,j+1} +\end{align} + +Assuming a constant $\alpha_x$ and $\alpha_y$ in each direction the equation can +be defined: + +#+NAME: eqn:genADI +\begin{align}\displaystyle +\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{n+1}_{ij}-C^{n+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n+1}_{ij}\right)}{\Delta y^2} +\end{align} + +Now we will define $s_x$ and $s_y$ respectively as followed: + +\begin{align}\displaystyle +s_x &= \frac{\alpha_x \cdot \frac{\Delta t}{2}}{\Delta x^2} \nonumber \\ +s_y &= \frac{\alpha_y \cdot \frac{\Delta t}{2}}{\Delta y^2} +\end{align} + +Equation [[eqn:genADI]], once developed in the $x$ direction, yields: + +\begin{align}\displaystyle +\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ +C^{n+1/2}_{ij}-C^n_{ij} &= s_x \delta^2_x C^{n+1/2}_{ij} + s_x \delta^2_y C^{n}_{ij} \nonumber \\ +-C^n_{ij} - s_x \delta^2_y C^{n}_{ij} &= C^{n+1/2}_{ij} + s_x \delta^2_x C^{n+1/2}_{ij} \nonumber \\ +-C^n_{ij} - s_x \left(C^n_{i,j-1}-2C^n_{i,j}+C^n_{i,j+1} \right)&= C^{n+1/2}_{ij} + s_x \left( C^{n+1/2}_{i-1,j} - 2C^{n+1/2}_{i,j} + C^{n+1/2}_{i+1,j} \right) +\end{align} + +All values on the left side of the equation are known and we are solving for +$C^{1/2}$. The calculation in $y$ direction for another $\frac{\Delta t}{2}$ is +similar and can also be easily derived. Therefore we do not show it here. + +*TODO*: +- ADI at boundaries #+LATEX: \clearpage From 7fd44aa67a64ac077240bca9b84edd4f549290cb Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 26 Apr 2022 10:07:51 +0200 Subject: [PATCH 16/36] Added 2D ADI scheme for boundaries --- doc/ADI_scheme.org | 39 +++++++++++++++++++++++++++------------ 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 3f3e303..431e0e5 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -186,12 +186,12 @@ implicit method. Therefore we make use of second difference operator defined in equation [[eqn:secondderivative]] in both $x$ and $y$ direction for half the time step $\Delta t$. We are denoting the numerator of equation [[eqn:1DBTCS]] as $\delta^2_x C^n_{ij}$ and $\delta^2_y C^n_{ij}$ respectively with $i$ the -position in $x$, $j$ the position in $y$ direction and $n$ as the current time +position in $x$, $j$ the position in $y$ direction and $T$ as the current time step: \begin{align}\displaystyle -\delta^2_x C^n_{ij} &= C^{n}_{i-1,j} - 2C^{n}_{i,j} + C^{n}_{i+1,j} \nonumber \\ -\delta^2_y C^n_{ij} &= C^{n}_{i,j-1} - 2C^{n}_{i,j} + C^{n}_{i,j+1} +\delta^2_x C^T_{ij} &= C^{T}_{i-1,j} - 2C^{T}_{i,j} + C^{T}_{i+1,j} \nonumber \\ +\delta^2_y C^T_{ij} &= C^{T}_{i,j-1} - 2C^{T}_{i,j} + C^{T}_{i,j+1} \end{align} Assuming a constant $\alpha_x$ and $\alpha_y$ in each direction the equation can @@ -199,8 +199,8 @@ be defined: #+NAME: eqn:genADI \begin{align}\displaystyle -\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ -\frac{C^{n+1}_{ij}-C^{n+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n+1}_{ij}\right)}{\Delta y^2} +\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{T+1}_{ij}-C^{T+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T+1}_{ij}\right)}{\Delta y^2} \end{align} Now we will define $s_x$ and $s_y$ respectively as followed: @@ -213,19 +213,34 @@ s_y &= \frac{\alpha_y \cdot \frac{\Delta t}{2}}{\Delta y^2} Equation [[eqn:genADI]], once developed in the $x$ direction, yields: \begin{align}\displaystyle -\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ -\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ -C^{n+1/2}_{ij}-C^n_{ij} &= s_x \delta^2_x C^{n+1/2}_{ij} + s_x \delta^2_y C^{n}_{ij} \nonumber \\ --C^n_{ij} - s_x \delta^2_y C^{n}_{ij} &= C^{n+1/2}_{ij} + s_x \delta^2_x C^{n+1/2}_{ij} \nonumber \\ --C^n_{ij} - s_x \left(C^n_{i,j-1}-2C^n_{i,j}+C^n_{i,j+1} \right)&= C^{n+1/2}_{ij} + s_x \left( C^{n+1/2}_{i-1,j} - 2C^{n+1/2}_{i,j} + C^{n+1/2}_{i+1,j} \right) +\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ +C^{T+1/2}_{ij}-C^T_{ij} &= s_x \delta^2_x C^{T+1/2}_{ij} + s_x \delta^2_y C^{T}_{ij} \nonumber \\ +-C^T_{ij} - s_x \delta^2_y C^{T}_{ij} &= C^{T+1/2}_{ij} + s_x \delta^2_x C^{T+1/2}_{ij} \end{align} All values on the left side of the equation are known and we are solving for $C^{1/2}$. The calculation in $y$ direction for another $\frac{\Delta t}{2}$ is similar and can also be easily derived. Therefore we do not show it here. +Substituting $\delta$ in the equation leads to following epxression: -*TODO*: -- ADI at boundaries +\begin{equation} + -C^T_{ij} - s_x \left(C^T_{i,j-1}-2C^T_{i,j}+C^T_{i,j+1} \right) = C^{T+1/2}_{ij} + s_x \left( C^{T+1/2}_{i-1,j} - 2C^{T+1/2}_{i,j} + C^{T+1/2}_{i+1,j} \right) +\end{equation} + +This scheme also only applies to inlet cells without a relation to boundaries. +Fortunately we already derived both cases of outer left and right inlet cell +respectively. Hence we are able to redefine each $\delta^2$ case in x and y +direction, assuming $l_x$ and $l_y$ the be the left boundary value and $r_x$ and +$r_y$ the right one for each direction $x$ and $y$. The equations are exemplary +for timestep $T+1/2$: + +\begin{align}\displaystyle +\delta^2_d C^{T+1/2}_{0,j} &= 2l_x - 3C^{T+1/2}_{0,j} + C^{T+1/2}_{1,j} \nonumber \\ +\delta^2_d C^{T+1/2}_{n,j} &= 2r_x - 3C^{T+1/2}_{n,j} + C^{T+1/2}_{n-1,j} \nonumber \\ +\delta^2_d C^{T+1/2}_{i,0} &= 2l_y - 3C^{T+1/2}_{i,0} + C^{T+1/2}_{i,1} \nonumber \\ +\delta^2_d C^{T+1/2}_{i,n} &= 2r_y - 3C^{T+1/2}_{i,n} + C^{T+1/2}_{i,n-1} +\end{align} #+LATEX: \clearpage From f42af515fb327117cdabe75a2df1612e52566508 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 10:06:38 +0200 Subject: [PATCH 17/36] Update applicable equation --- doc/ADI_scheme.org | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 8a3fd92..557fe34 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -163,11 +163,17 @@ C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) \end{align} -Which leads to the following term applicable to our approach: +Now we define variable $s_x$ as following: + +\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: #+NAME: eqn:5 \begin{equation}\displaystyle - -C^j_0} = \left[\frac{\alpha \cdot \Delta t}{\Delta x^2}\right] \cdot \left(C^{j+1}_1 + 2l \right) + \left[ -1 - 3 \cdot \frac{\alpha \cdot \Delta t}{\Delta x^2} \right] \cdot C^{j+1}_0 + -C^j_0} = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} In case of constant right boundary, the finite difference of point From f31937f4e0cff661611cb318b829e8d99f3a08b6 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 10:09:43 +0200 Subject: [PATCH 18/36] Fix broken reference --- doc/ADI_scheme.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 557fe34..3cabbf1 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -192,7 +192,7 @@ C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{n-1} \right) \end{align} -If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:5]] +If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:6]] becomes zero and we are left with: From 29d3d81ff4e41e649fe05ec81402d5a848326c36 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 11:27:30 +0200 Subject: [PATCH 19/36] Apply bc to ghost node --- app/main_1D.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/app/main_1D.cpp b/app/main_1D.cpp index e908039..2fb7c21 100644 --- a/app/main_1D.cpp +++ b/app/main_1D.cpp @@ -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[1] = {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)); From 08da0b47f44159f208319b81d9b995e79eed25e9 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 11:30:59 +0200 Subject: [PATCH 20/36] Revert to 'age' state and append with implicit BTCS --- doc/ADI_scheme.org | 96 ++++++++++++++++++++++++++++++---------------- 1 file changed, 64 insertions(+), 32 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 3cabbf1..09158ed 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -132,7 +132,7 @@ and assuming constant $\alpha$: #+NAME: eqn:2 \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} + \frac{C_i^{j+1} -C_i^{j}}{\Delta t} = \alpha\frac{\frac{C^j_{i+1}-C^j_{i}}{\Delta x}-\frac{C^j_{i}-C^j_{i-1}}{\Delta x}}{\Delta x} \end{equation} In practice, we evaluate the first derivatives of $C$ w.r.t. $x$ on @@ -150,6 +150,69 @@ 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^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^j_{1}-C^j_{0}}{\Delta x}- +\frac{C^j_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +#+NAME: eqn:4 +\begin{align}\displaystyle +C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^j_{1}-C^j_{0}- 2 C^j_{0}+2l \right) \nonumber \\ + & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^j_{1}- 3 C^j_{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^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}- +\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +Which, developed, gives +#+NAME: eqn:6 +\begin{align}\displaystyle +C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\ + & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{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^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n) +\end{equation} + + + +A similar treatment can be applied to the BTCS implicit scheme. + +*** implicit BTCS + +\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} + +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 +boundaries. 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: + \begin{equation}\displaystyle \frac{C_0^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}- \frac{C^{j+1}_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} @@ -157,7 +220,6 @@ calling $l$ the numerical value of a constant boundary condition: This expression, once developed, yields: -#+NAME: eqn:4 \begin{align}\displaystyle C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}-C^{j+1}_{0}- 2 C^{j+1}_{0}+2l \right) \nonumber \\ & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) @@ -171,36 +233,6 @@ Now we define variable $s_x$ as following: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: -#+NAME: eqn:5 \begin{equation}\displaystyle -C^j_0} = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} - -In case of constant right boundary, the finite difference of point -$C_n$ - calling $r$ the right boundary value - is: - -#+NAME: eqn:6 -\begin{equation}\displaystyle -\frac{C_n^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}- -\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x} -\end{equation} - -Which, developed, gives -#+NAME: eqn:7 -\begin{align}\displaystyle -C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\ - & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{n-1} \right) -\end{align} - -If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:6]] -becomes zero and we are left with: - - -#+NAME: eqn:8 -\begin{equation}\displaystyle -C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n) -\end{equation} - - - -A similar treatment can be applied to the BTCS implicit scheme. From 0f646d478c63a780817eaded175a4d2a361677c9 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Thu, 21 Apr 2022 12:51:34 +0200 Subject: [PATCH 21/36] Fix unnecessary bracket --- doc/ADI_scheme.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 09158ed..d47f3cb 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -234,5 +234,5 @@ Now we define variable $s_x$ as following: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: \begin{equation}\displaystyle - -C^j_0} = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + -C^j_0 = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} From d54fe25cace10879e7846f67cedc20b3d8955db1 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Fri, 22 Apr 2022 09:49:53 +0200 Subject: [PATCH 22/36] Update documentation of implicit BTCS --- doc/ADI_scheme.org | 41 +++++++++++++++++++++++++++-------------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index d47f3cb..7a0b0c7 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -195,23 +195,36 @@ A similar treatment can be applied to the BTCS implicit scheme. *** implicit BTCS -\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} +First, we define the Backward time difference: + +\begin{equation} + \frac{\partial C }{\partial t} = \frac{C^j_i - C^{j-1}_i}{\Delta t} \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$. +Second the spatial derivative approximation: -This discretization works for all internal cells, but not for the -boundaries. To properly treat them, we need to account for the -discrepancy in the discretization. +\begin{equation} + \frac{\partial^2 C }{\partial t} = \frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{i-1}}{\Delta x}}{\Delta x} +\end{equation} -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: +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} -C_i^{j-1}}{\Delta t} = \alpha\frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{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: + +\begin{align}\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} \nonumber \\ + & = \alpha\frac{C^{j+1}_{i-1} - 2C^{j+1}_{i} + C^{j+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^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}- @@ -225,7 +238,7 @@ C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) \end{align} -Now we define variable $s_x$ as following: +Now we define variable $s_x$ as followed: \begin{equation} s_x = \frac{\alpha \cdot \Delta t}{\Delta x^2} From 8b4f1aae46164c707f52c1cc7e2c6732467c78bd Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Fri, 22 Apr 2022 13:55:50 +0200 Subject: [PATCH 23/36] MDL: some restructuring in ADI_scheme.org --- doc/ADI_scheme.org | 259 ++++++++++++++++++++++++--------------------- 1 file changed, 136 insertions(+), 123 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 7a0b0c7..a3497aa 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -1,114 +1,11 @@ -#+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} #+OPTIONS: toc:nil -* Input -- =c= $\rightarrow c$ - - containing current concentrations at each grid cell for species - - size: $N \times M$ - - row-major -- =alpha= $\rightarrow \alpha$ - - diffusion coefficient for both directions (x and y) - - size: $N \times M$ - - row-major -- =boundary_condition= $\rightarrow bc$ - - Defines closed or constant boundary condition for each grid cell - - size: $N \times M$ - - row-major - -* Internals - -- =A_matrix= $\rightarrow A$ - - coefficient matrix for linear equation system implemented as sparse matrix - - size: $((N+2)\cdot M) \times ((N+2)\cdot M)$ (including ghost zones in x direction) - - column-major (not relevant) - -- =b_vector= $\rightarrow b$ - - right hand side of the linear equation system - - size: $(N+2) \cdot M$ - - column-major (not relevant) -- =x_vector= $\rightarrow x$ - - solutions of the linear equation system - - size: $(N+2) \cdot M$ - - column-major (not relevant) - -* Calculation for $\frac{1}{2}$ timestep - -** Symbolic addressing of grid cells -[[./grid.png]] - -** Filling of matrix $A$ - -- row-wise iterating with $i$ over =c= and =\alpha= matrix respectively -- addressing each element of a row with $j$ -- matrix $A$ also containing $+2$ ghost nodes for each row of input matrix $\alpha$ - - $\rightarrow offset = N+2$ - - addressing each object $(i,j)$ in matrix $A$ with $(offset \cdot i + j, offset \cdot i + j)$ - -*** Rules - -$s_x(i,j) = \frac{\alpha(i,j)*\frac{t}{2}}{\Delta x^2}$ where $x$ defining the domain size in x direction. - -For the sake of simplicity we assume that each row of the $A$ matrix is addressed correctly with the given offset. - -**** Ghost nodes - -$A(i,-1) = 1$ - -$A(i,N) = 1$ - -**** Inlet - -$A(i,j) = \begin{cases} -1 & \text{if } bc(i,j) = \text{constant} \\ --1-2*s_x(i,j) & \text{else} -\end{cases}$ - -$A(i,j\pm 1) = \begin{cases} -0 & \text{if } bc(i,j) = \text{constant} \\ -s_x(i,j) & \text{else} -\end{cases}$ - -** Filling of vector $b$ - -- each elements assign a concrete value to the according value of the row of matrix $A$ -- 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 - -$b(i,-1) = \begin{cases} -0 & \text{if } bc(i,0) = \text{constant} \\ -c(i,0) & \text{else} -\end{cases}$ - -$b(i,N) = \begin{cases} -0 & \text{if } bc(i,N-1) = \text{constant} \\ -c(i,N-1) & \text{else} -\end{cases}$ - -*** 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] - -$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 - -** Finite differences with nodes as cells' centres - -*** The explicit FTCS scheme as in PHREEQC +* Finite differences with nodes as cells' centres The 1D diffusion equation is: @@ -118,17 +15,22 @@ The 1D diffusion equation is: & = \alpha \frac{\partial^2 C}{\partial x^2} \end{align} -We discretize it following a Forward Time, Centered Space finite -difference scheme where the nodes correspond to the centers of a grid -such as: +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 is in $x=dx/2$, with $dx=L/n$. +cell - which are the points constituting the finite difference nodes - +is in $x=dx/2$, with $dx=L/n$. -We discretize [[eqn:1]] as following, for each index i in 1, \dots, n-1 -and assuming constant $\alpha$: + +** 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 @@ -142,8 +44,8 @@ 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 -boundaries. To properly treat them, we need to account for the -discrepancy in the discretization. +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, @@ -193,29 +95,30 @@ C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} A similar treatment can be applied to the BTCS implicit scheme. -*** implicit BTCS +** Implicit BTCS scheme -First, we define the Backward time difference: +First, we define the Backward Time difference: \begin{equation} - \frac{\partial C }{\partial t} = \frac{C^j_i - C^{j-1}_i}{\Delta t} + \frac{\partial C^{j+1} }{\partial t} = \frac{C^{j+1}_i - C^{j}_i}{\Delta t} \end{equation} -Second the spatial derivative approximation: +Second the spatial derivative approximation, evaluated at time level $j+1$: \begin{equation} - \frac{\partial^2 C }{\partial t} = \frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{i-1}}{\Delta x}}{\Delta x} + \frac{\partial^2 C^{j+1} }{\partial x^2} = \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} 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} -C_i^{j-1}}{\Delta t} = \alpha\frac{\frac{C^{j}_{i+1}-C^{j}_{i}}{\Delta x}-\frac{C^{j}_{i}-C^{j}_{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: +# \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: \begin{align}\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} \nonumber \\ @@ -249,3 +152,113 @@ Substituting with the new variable $s_x$ and reordering of terms leads to the eq \begin{equation}\displaystyle -C^j_0 = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 \end{equation} + +*TODO* +- Right boundary +- Tridiagonal matrix filling + + + + +#+LATEX: \clearpage + +* Old stuff + +** Input + +- =c= $\rightarrow c$ + - containing current concentrations at each grid cell for species + - size: $N \times M$ + - row-major +- =alpha= $\rightarrow \alpha$ + - diffusion coefficient for both directions (x and y) + - size: $N \times M$ + - row-major +- =boundary_condition= $\rightarrow bc$ + - Defines closed or constant boundary condition for each grid cell + - size: $N \times M$ + - row-major + +** Internals + +- =A_matrix= $\rightarrow A$ + - coefficient matrix for linear equation system implemented as sparse matrix + - size: $((N+2)\cdot M) \times ((N+2)\cdot M)$ (including ghost zones in x direction) + - column-major (not relevant) + +- =b_vector= $\rightarrow b$ + - right hand side of the linear equation system + - size: $(N+2) \cdot M$ + - column-major (not relevant) +- =x_vector= $\rightarrow x$ + - solutions of the linear equation system + - size: $(N+2) \cdot M$ + - column-major (not relevant) + +** Calculation for $\frac{1}{2}$ timestep + +** Symbolic addressing of grid cells +[[./grid.png]] + +** Filling of matrix $A$ + +- row-wise iterating with $i$ over =c= and =\alpha= matrix respectively +- addressing each element of a row with $j$ +- matrix $A$ also containing $+2$ ghost nodes for each row of input matrix $\alpha$ + - $\rightarrow offset = N+2$ + - addressing each object $(i,j)$ in matrix $A$ with $(offset \cdot i + j, offset \cdot i + j)$ + +*** Rules + +$s_x(i,j) = \frac{\alpha(i,j)*\frac{t}{2}}{\Delta x^2}$ where $x$ defining the domain size in x direction. + +For the sake of simplicity we assume that each row of the $A$ matrix is addressed correctly with the given offset. + +**** Ghost nodes + +$A(i,-1) = 1$ + +$A(i,N) = 1$ + +**** Inlet + +$A(i,j) = \begin{cases} +1 & \text{if } bc(i,j) = \text{constant} \\ +-1-2*s_x(i,j) & \text{else} +\end{cases}$ + +$A(i,j\pm 1) = \begin{cases} +0 & \text{if } bc(i,j) = \text{constant} \\ +s_x(i,j) & \text{else} +\end{cases}$ + +** Filling of vector $b$ + +- each elements assign a concrete value to the according value of the row of matrix $A$ +- 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 + +$b(i,-1) = \begin{cases} +0 & \text{if } bc(i,0) = \text{constant} \\ +c(i,0) & \text{else} +\end{cases}$ + +$b(i,N) = \begin{cases} +0 & \text{if } bc(i,N-1) = \text{constant} \\ +c(i,N-1) & \text{else} +\end{cases}$ + +*** 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}$ + +\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}$ From aa91824eefe82eff854240d584df6b58e38c8216 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 25 Apr 2022 10:52:23 +0200 Subject: [PATCH 24/36] Reorder equation for better readability --- doc/ADI_scheme.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index a3497aa..d363194 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -150,7 +150,7 @@ Now we define variable $s_x$ as followed: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: \begin{equation}\displaystyle - -C^j_0 = s_x \cdot C^{j+1}_1 + (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + -C^j_0 = (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + s_x \cdot C^{j+1}_1 \end{equation} *TODO* From e8b1c1daae4395dc29b36181cc8f1d8257a22ad7 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 25 Apr 2022 10:50:06 +0200 Subject: [PATCH 25/36] Added right boundary --- doc/ADI_scheme.org | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index d363194..391d3a6 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -153,8 +153,27 @@ Substituting with the new variable $s_x$ and reordering of terms leads to the eq -C^j_0 = (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + s_x \cdot C^{j+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^{j+1} -C_n^{j}}{\Delta t} = \alpha\frac{\frac{r-C^{j+1}_{n}}{\frac{\Delta x}{2}}- +\frac{C^{j+1}_{n}-C^{j+1}_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +\begin{align}\displaystyle +C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 2C^{j+1}_{n} - C^{j+1}_{n} + C^{j+1}_{n-1} \right) \nonumber \\ + & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 3C^{j+1}_{n} + C^{j+1}_{n-1} \right) +\end{align} + +Now rearrange terms and substituting with $s_x$ leads to: + +\begin{equation}\displaystyle + -C^j_n = s_x \cdot C^{j+1}_{n-1} + (-1 - 3s_x) \cdot C^{j+1}_n + (2s_x) \cdot r +\end{equation} + *TODO* -- Right boundary - Tridiagonal matrix filling From ffefb6ab50339dde222e27ae4a44b47cbcc9f0e8 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Mon, 25 Apr 2022 12:31:27 +0200 Subject: [PATCH 26/36] Added 2D ADI scheme --- doc/ADI_scheme.org | 48 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 391d3a6..3f3e303 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -105,6 +105,7 @@ First, we define the Backward Time difference: Second the spatial derivative approximation, evaluated at time level $j+1$: +#+NAME: eqn:secondderivative \begin{equation} \frac{\partial^2 C^{j+1} }{\partial x^2} = \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} @@ -120,6 +121,7 @@ equations given above leads to the following 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^{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} \nonumber \\ & = \alpha\frac{C^{j+1}_{i-1} - 2C^{j+1}_{i} + C^{j+1}_{i+1}}{\Delta x^2} @@ -176,8 +178,54 @@ Now rearrange terms and substituting with $s_x$ leads to: *TODO* - Tridiagonal matrix filling +** 2D ADI using BTCS scheme +In the previous sections we described the usage of FTCS and BTCS on 1D grids. To +make usage of the BTCS scheme we are following the alternating-direction +implicit method. Therefore we make use of second difference operator defined in +equation [[eqn:secondderivative]] in both $x$ and $y$ direction for half the time +step $\Delta t$. We are denoting the numerator of equation [[eqn:1DBTCS]] as +$\delta^2_x C^n_{ij}$ and $\delta^2_y C^n_{ij}$ respectively with $i$ the +position in $x$, $j$ the position in $y$ direction and $n$ as the current time +step: +\begin{align}\displaystyle +\delta^2_x C^n_{ij} &= C^{n}_{i-1,j} - 2C^{n}_{i,j} + C^{n}_{i+1,j} \nonumber \\ +\delta^2_y C^n_{ij} &= C^{n}_{i,j-1} - 2C^{n}_{i,j} + C^{n}_{i,j+1} +\end{align} + +Assuming a constant $\alpha_x$ and $\alpha_y$ in each direction the equation can +be defined: + +#+NAME: eqn:genADI +\begin{align}\displaystyle +\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{n+1}_{ij}-C^{n+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n+1}_{ij}\right)}{\Delta y^2} +\end{align} + +Now we will define $s_x$ and $s_y$ respectively as followed: + +\begin{align}\displaystyle +s_x &= \frac{\alpha_x \cdot \frac{\Delta t}{2}}{\Delta x^2} \nonumber \\ +s_y &= \frac{\alpha_y \cdot \frac{\Delta t}{2}}{\Delta y^2} +\end{align} + +Equation [[eqn:genADI]], once developed in the $x$ direction, yields: + +\begin{align}\displaystyle +\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ +C^{n+1/2}_{ij}-C^n_{ij} &= s_x \delta^2_x C^{n+1/2}_{ij} + s_x \delta^2_y C^{n}_{ij} \nonumber \\ +-C^n_{ij} - s_x \delta^2_y C^{n}_{ij} &= C^{n+1/2}_{ij} + s_x \delta^2_x C^{n+1/2}_{ij} \nonumber \\ +-C^n_{ij} - s_x \left(C^n_{i,j-1}-2C^n_{i,j}+C^n_{i,j+1} \right)&= C^{n+1/2}_{ij} + s_x \left( C^{n+1/2}_{i-1,j} - 2C^{n+1/2}_{i,j} + C^{n+1/2}_{i+1,j} \right) +\end{align} + +All values on the left side of the equation are known and we are solving for +$C^{1/2}$. The calculation in $y$ direction for another $\frac{\Delta t}{2}$ is +similar and can also be easily derived. Therefore we do not show it here. + +*TODO*: +- ADI at boundaries #+LATEX: \clearpage From 57cfa682c576311a3404b1bb3e127811cc7bc722 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Tue, 26 Apr 2022 10:07:51 +0200 Subject: [PATCH 27/36] Added 2D ADI scheme for boundaries --- doc/ADI_scheme.org | 39 +++++++++++++++++++++++++++------------ 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 3f3e303..431e0e5 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -186,12 +186,12 @@ implicit method. Therefore we make use of second difference operator defined in equation [[eqn:secondderivative]] in both $x$ and $y$ direction for half the time step $\Delta t$. We are denoting the numerator of equation [[eqn:1DBTCS]] as $\delta^2_x C^n_{ij}$ and $\delta^2_y C^n_{ij}$ respectively with $i$ the -position in $x$, $j$ the position in $y$ direction and $n$ as the current time +position in $x$, $j$ the position in $y$ direction and $T$ as the current time step: \begin{align}\displaystyle -\delta^2_x C^n_{ij} &= C^{n}_{i-1,j} - 2C^{n}_{i,j} + C^{n}_{i+1,j} \nonumber \\ -\delta^2_y C^n_{ij} &= C^{n}_{i,j-1} - 2C^{n}_{i,j} + C^{n}_{i,j+1} +\delta^2_x C^T_{ij} &= C^{T}_{i-1,j} - 2C^{T}_{i,j} + C^{T}_{i+1,j} \nonumber \\ +\delta^2_y C^T_{ij} &= C^{T}_{i,j-1} - 2C^{T}_{i,j} + C^{T}_{i,j+1} \end{align} Assuming a constant $\alpha_x$ and $\alpha_y$ in each direction the equation can @@ -199,8 +199,8 @@ be defined: #+NAME: eqn:genADI \begin{align}\displaystyle -\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ -\frac{C^{n+1}_{ij}-C^{n+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n+1}_{ij}\right)}{\Delta y^2} +\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{T+1}_{ij}-C^{T+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T+1}_{ij}\right)}{\Delta y^2} \end{align} Now we will define $s_x$ and $s_y$ respectively as followed: @@ -213,19 +213,34 @@ s_y &= \frac{\alpha_y \cdot \frac{\Delta t}{2}}{\Delta y^2} Equation [[eqn:genADI]], once developed in the $x$ direction, yields: \begin{align}\displaystyle -\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} + \delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ -\frac{C^{n+1/2}_{ij}-C^n_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{n+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{n}_{ij}\right)}{\Delta x^2} \nonumber \\ -C^{n+1/2}_{ij}-C^n_{ij} &= s_x \delta^2_x C^{n+1/2}_{ij} + s_x \delta^2_y C^{n}_{ij} \nonumber \\ --C^n_{ij} - s_x \delta^2_y C^{n}_{ij} &= C^{n+1/2}_{ij} + s_x \delta^2_x C^{n+1/2}_{ij} \nonumber \\ --C^n_{ij} - s_x \left(C^n_{i,j-1}-2C^n_{i,j}+C^n_{i,j+1} \right)&= C^{n+1/2}_{ij} + s_x \left( C^{n+1/2}_{i-1,j} - 2C^{n+1/2}_{i,j} + C^{n+1/2}_{i+1,j} \right) +\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ +\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ +C^{T+1/2}_{ij}-C^T_{ij} &= s_x \delta^2_x C^{T+1/2}_{ij} + s_x \delta^2_y C^{T}_{ij} \nonumber \\ +-C^T_{ij} - s_x \delta^2_y C^{T}_{ij} &= C^{T+1/2}_{ij} + s_x \delta^2_x C^{T+1/2}_{ij} \end{align} All values on the left side of the equation are known and we are solving for $C^{1/2}$. The calculation in $y$ direction for another $\frac{\Delta t}{2}$ is similar and can also be easily derived. Therefore we do not show it here. +Substituting $\delta$ in the equation leads to following epxression: -*TODO*: -- ADI at boundaries +\begin{equation} + -C^T_{ij} - s_x \left(C^T_{i,j-1}-2C^T_{i,j}+C^T_{i,j+1} \right) = C^{T+1/2}_{ij} + s_x \left( C^{T+1/2}_{i-1,j} - 2C^{T+1/2}_{i,j} + C^{T+1/2}_{i+1,j} \right) +\end{equation} + +This scheme also only applies to inlet cells without a relation to boundaries. +Fortunately we already derived both cases of outer left and right inlet cell +respectively. Hence we are able to redefine each $\delta^2$ case in x and y +direction, assuming $l_x$ and $l_y$ the be the left boundary value and $r_x$ and +$r_y$ the right one for each direction $x$ and $y$. The equations are exemplary +for timestep $T+1/2$: + +\begin{align}\displaystyle +\delta^2_d C^{T+1/2}_{0,j} &= 2l_x - 3C^{T+1/2}_{0,j} + C^{T+1/2}_{1,j} \nonumber \\ +\delta^2_d C^{T+1/2}_{n,j} &= 2r_x - 3C^{T+1/2}_{n,j} + C^{T+1/2}_{n-1,j} \nonumber \\ +\delta^2_d C^{T+1/2}_{i,0} &= 2l_y - 3C^{T+1/2}_{i,0} + C^{T+1/2}_{i,1} \nonumber \\ +\delta^2_d C^{T+1/2}_{i,n} &= 2r_y - 3C^{T+1/2}_{i,n} + C^{T+1/2}_{i,n-1} +\end{align} #+LATEX: \clearpage From e9a1d067849a2b5dcb49512458c3c4bc005121d6 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 09:55:38 +0200 Subject: [PATCH 28/36] Apply new scheme to model (only 1D) --- src/BTCSDiffusion.cpp | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 1271156..606d6eb 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -239,8 +239,8 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( A_matrix.insert(1, 1) = 1; } else { double sx = (alpha[0] * time_step) / (dx * dx); - A_matrix.insert(1, 1) = -1. - 2. * sx; - A_matrix.insert(1, 0) = sx; + A_matrix.insert(1, 1) = -1. - 3. * sx; + A_matrix.insert(1, 0) = 2. * sx; A_matrix.insert(1, 2) = sx; } @@ -261,9 +261,9 @@ void Diffusion::BTCSDiffusion::fillMatrixFromRow( 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. - 2. * sx; + 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) = sx; + A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx; } A_matrix.insert(A_size - 1, A_size - 1) = 1; @@ -294,15 +294,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) { From d35a27f54a4adfa6b8b1ca5e59302f24b51ff9d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 26 Apr 2022 13:29:12 +0200 Subject: [PATCH 29/36] Apply 2D scheme to model --- src/BTCSDiffusion.cpp | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 606d6eb..4f7b109 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -141,11 +141,10 @@ 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++) { @@ -184,12 +183,17 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c, // 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 @@ -210,12 +214,17 @@ 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; From 0211e3d9b2adf36969f0281ae882abe0d05fd0c6 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Tue, 26 Apr 2022 17:47:57 +0200 Subject: [PATCH 30/36] Improved consistency of notation in ADI_scheme and reworked section 2D ADI --- doc/ADI_scheme.org | 188 ++++++++++++++++++++++++++------------------- 1 file changed, 108 insertions(+), 80 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 431e0e5..dac13bd 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -1,11 +1,13 @@ #+TITLE: Numerical solution of diffusion equation in 2D with ADI Scheme #+LaTeX_CLASS_OPTIONS: [a4paper,10pt] #+LATEX_HEADER: \usepackage{fullpage} -#+LATEX_HEADER: \usepackage{amsmath} +#+LATEX_HEADER: \usepackage{amsmath, systeme} #+OPTIONS: toc:nil -* Finite differences with nodes as cells' centres +* Diffusion in 1D + +** Finite differences with nodes as cells' centres The 1D diffusion equation is: @@ -34,7 +36,7 @@ $\alpha$, we can write: #+NAME: eqn:2 \begin{equation}\displaystyle - \frac{C_i^{j+1} -C_i^{j}}{\Delta t} = \alpha\frac{\frac{C^j_{i+1}-C^j_{i}}{\Delta x}-\frac{C^j_{i}-C^j_{i-1}}{\Delta x}}{\Delta x} + \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 @@ -53,16 +55,16 @@ calling $l$ the numerical value of a constant boundary condition: #+NAME: eqn:3 \begin{equation}\displaystyle -\frac{C_0^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^j_{1}-C^j_{0}}{\Delta x}- -\frac{C^j_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} +\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^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^j_{1}-C^j_{0}- 2 C^j_{0}+2l \right) \nonumber \\ - & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^j_{1}- 3 C^j_{0} +2l \right) +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} @@ -71,15 +73,15 @@ $C_n$ - calling $r$ the right boundary value - is: #+NAME: eqn:5 \begin{equation}\displaystyle -\frac{C_n^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}- -\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x} +\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^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\ - & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{n-1} \right) +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]] @@ -88,7 +90,7 @@ becomes zero and we are left with: #+NAME: eqn:7 \begin{equation}\displaystyle -C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n) +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} @@ -100,14 +102,14 @@ A similar treatment can be applied to the BTCS implicit scheme. First, we define the Backward Time difference: \begin{equation} - \frac{\partial C^{j+1} }{\partial t} = \frac{C^{j+1}_i - C^{j}_i}{\Delta t} + \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 $j+1$: +Second the spatial derivative approximation, evaluated at time level $t+1$: #+NAME: eqn:secondderivative \begin{equation} - \frac{\partial^2 C^{j+1} }{\partial x^2} = \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} + \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 @@ -123,8 +125,8 @@ equations given above leads to the following equation: #+NAME: eqn:1DBTCS \begin{align}\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} \nonumber \\ - & = \alpha\frac{C^{j+1}_{i-1} - 2C^{j+1}_{i} + C^{j+1}_{i+1}}{\Delta x^2} +\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 @@ -132,15 +134,15 @@ 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^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}- -\frac{C^{j+1}_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} +\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^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}-C^{j+1}_{0}- 2 C^{j+1}_{0}+2l \right) \nonumber \\ - & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) +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: @@ -152,96 +154,122 @@ Now we define variable $s_x$ as followed: Substituting with the new variable $s_x$ and reordering of terms leads to the equation applicable to our model: \begin{equation}\displaystyle - -C^j_0 = (2s_x) \cdot l + (-1 - 3s_x) \cdot C^{j+1}_0 + s_x \cdot C^{j+1}_1 + -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^{j+1} -C_n^{j}}{\Delta t} = \alpha\frac{\frac{r-C^{j+1}_{n}}{\frac{\Delta x}{2}}- -\frac{C^{j+1}_{n}-C^{j+1}_{n-1}}{\Delta x}}{\Delta x} +\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^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 2C^{j+1}_{n} - C^{j+1}_{n} + C^{j+1}_{n-1} \right) \nonumber \\ - & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2r - 3C^{j+1}_{n} + C^{j+1}_{n-1} \right) +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^j_n = s_x \cdot C^{j+1}_{n-1} + (-1 - 3s_x) \cdot C^{j+1}_n + (2s_x) \cdot r + -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 -** 2D ADI using BTCS scheme +#+LATEX: \clearpage -In the previous sections we described the usage of FTCS and BTCS on 1D grids. To -make usage of the BTCS scheme we are following the alternating-direction -implicit method. Therefore we make use of second difference operator defined in -equation [[eqn:secondderivative]] in both $x$ and $y$ direction for half the time -step $\Delta t$. We are denoting the numerator of equation [[eqn:1DBTCS]] as -$\delta^2_x C^n_{ij}$ and $\delta^2_y C^n_{ij}$ respectively with $i$ the -position in $x$, $j$ the position in $y$ direction and $T$ as the current time -step: +* Diffusion in 2D: the Alternating Direction Implicit scheme -\begin{align}\displaystyle -\delta^2_x C^T_{ij} &= C^{T}_{i-1,j} - 2C^{T}_{i,j} + C^{T}_{i+1,j} \nonumber \\ -\delta^2_y C^T_{ij} &= C^{T}_{i,j-1} - 2C^{T}_{i,j} + C^{T}_{i,j+1} -\end{align} -Assuming a constant $\alpha_x$ and $\alpha_y$ in each direction the equation can -be defined: - -#+NAME: eqn:genADI -\begin{align}\displaystyle -\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ -\frac{C^{T+1}_{ij}-C^{T+1/2}_{ij}}{\frac{\Delta t}{2}} &= \alpha_y \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T+1}_{ij}\right)}{\Delta y^2} -\end{align} - -Now we will define $s_x$ and $s_y$ respectively as followed: - -\begin{align}\displaystyle -s_x &= \frac{\alpha_x \cdot \frac{\Delta t}{2}}{\Delta x^2} \nonumber \\ -s_y &= \frac{\alpha_y \cdot \frac{\Delta t}{2}}{\Delta y^2} -\end{align} - -Equation [[eqn:genADI]], once developed in the $x$ direction, yields: - -\begin{align}\displaystyle -\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} + \delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ -\frac{C^{T+1/2}_{ij}-C^T_{ij}}{\frac{\Delta t}{2}} &= \alpha_x \frac{\left( \delta^2_x C^{T+1/2}_{ij} \right)}{\Delta x^2} + \alpha_x \frac{\left(\delta^2_y C^{T}_{ij}\right)}{\Delta x^2} \nonumber \\ -C^{T+1/2}_{ij}-C^T_{ij} &= s_x \delta^2_x C^{T+1/2}_{ij} + s_x \delta^2_y C^{T}_{ij} \nonumber \\ --C^T_{ij} - s_x \delta^2_y C^{T}_{ij} &= C^{T+1/2}_{ij} + s_x \delta^2_x C^{T+1/2}_{ij} -\end{align} - -All values on the left side of the equation are known and we are solving for -$C^{1/2}$. The calculation in $y$ direction for another $\frac{\Delta t}{2}$ is -similar and can also be easily derived. Therefore we do not show it here. -Substituting $\delta$ in the equation leads to following epxression: +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} - -C^T_{ij} - s_x \left(C^T_{i,j-1}-2C^T_{i,j}+C^T_{i,j+1} \right) = C^{T+1/2}_{ij} + s_x \left( C^{T+1/2}_{i-1,j} - 2C^{T+1/2}_{i,j} + C^{T+1/2}_{i+1,j} \right) +\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} -This scheme also only applies to inlet cells without a relation to boundaries. -Fortunately we already derived both cases of outer left and right inlet cell -respectively. Hence we are able to redefine each $\delta^2$ case in x and y -direction, assuming $l_x$ and $l_y$ the be the left boundary value and $r_x$ and -$r_y$ the right one for each direction $x$ and $y$. The equations are exemplary -for timestep $T+1/2$: +** 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 -\delta^2_d C^{T+1/2}_{0,j} &= 2l_x - 3C^{T+1/2}_{0,j} + C^{T+1/2}_{1,j} \nonumber \\ -\delta^2_d C^{T+1/2}_{n,j} &= 2r_x - 3C^{T+1/2}_{n,j} + C^{T+1/2}_{n-1,j} \nonumber \\ -\delta^2_d C^{T+1/2}_{i,0} &= 2l_y - 3C^{T+1/2}_{i,0} + C^{T+1/2}_{i,1} \nonumber \\ -\delta^2_d C^{T+1/2}_{i,n} &= 2r_y - 3C^{T+1/2}_{i,n} + C^{T+1/2}_{i,n-1} +\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 inlet cells without a relation to +boundaries. Fortunately we already derived both cases of outer left +and right inlet cell respectively. Hence we are able to redefine each +$\delta^2$ case in x and y direction, assuming $l_x$ and $l_y$ the be +the left boundary value and $r_x$ and $r_y$ the right one for each +direction $x$ and $y$. The equations are exemplary for time level +$t+1/2$: + +\begin{equation} +\systeme{ + \displaystyle \delta^2_d C^{t+1/2}_{0,j} = 2l_x - 3C^{t+1/2}_{0,j} + C^{t+1/2}_{1,j} , + \displaystyle \delta^2_d C^{t+1/2}_{n,j} = 2r_x - 3C^{t+1/2}_{n,j} + C^{t+1/2}_{n-1,j} , + \displaystyle \delta^2_d C^{t+1/2}_{i,0} = 2l_y - 3C^{t+1/2}_{i,0} + C^{t+1/2}_{i,1} , + \displaystyle \delta^2_d C^{t+1/2}_{i,n} = 2r_y - 3C^{t+1/2}_{i,n} + C^{t+1/2}_{i,n-1} +} +\end{equation} + #+LATEX: \clearpage * Old stuff From 45f080271514678bb8d28dd3baaab8182afda0e3 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Tue, 26 Apr 2022 19:08:32 +0200 Subject: [PATCH 31/36] Improved docs about boundary conditions treatment for 2D ADI --- doc/ADI_scheme.org | 56 +++++++++++++++++++++++++++++++++++++--------- 1 file changed, 45 insertions(+), 11 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index dac13bd..50544fe 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -253,23 +253,57 @@ yields: -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 inlet cells without a relation to -boundaries. Fortunately we already derived both cases of outer left -and right inlet cell respectively. Hence we are able to redefine each -$\delta^2$ case in x and y direction, assuming $l_x$ and $l_y$ the be -the left boundary value and $r_x$ and $r_y$ the right one for each -direction $x$ and $y$. The equations are exemplary for time level -$t+1/2$: +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 \delta^2_d C^{t+1/2}_{0,j} = 2l_x - 3C^{t+1/2}_{0,j} + C^{t+1/2}_{1,j} , - \displaystyle \delta^2_d C^{t+1/2}_{n,j} = 2r_x - 3C^{t+1/2}_{n,j} + C^{t+1/2}_{n-1,j} , - \displaystyle \delta^2_d C^{t+1/2}_{i,0} = 2l_y - 3C^{t+1/2}_{i,0} + C^{t+1/2}_{i,1} , - \displaystyle \delta^2_d C^{t+1/2}_{i,n} = 2r_y - 3C^{t+1/2}_{i,n} + C^{t+1/2}_{i,n-1} + \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 From 628e3288beac165ef72a80c27c7e8a67619676c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 27 Apr 2022 12:37:43 +0200 Subject: [PATCH 32/36] Added 2D example with constant bc on left side --- app/CMakeLists.txt | 3 +++ app/main_2D_const.cpp | 57 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+) create mode 100644 app/main_2D_const.cpp 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_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; +} From 9706d9a4b170e631b6adcb959d457cd37c0ff59b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 27 Apr 2022 12:38:40 +0200 Subject: [PATCH 33/36] Fix indexiation of bc field --- src/BTCSDiffusion.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp index 4f7b109..54b9789 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -149,8 +149,8 @@ void Diffusion::BTCSDiffusion::simulate2D( #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; } @@ -162,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(); } } @@ -183,7 +183,7 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c, // first, iterate over first row for (int j = 0; j < n_cols; j++) { - boundary_condition tmp_bc = bc(0,j+1); + boundary_condition tmp_bc = bc(0, j + 1); if (tmp_bc.type == Diffusion::BC_CLOSED) continue; @@ -193,7 +193,7 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c, 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); + (2 * y_values[0] - 3 * y_values[1] + y_values[2]) / (dx * dx); } // then iterate over inlet @@ -214,7 +214,7 @@ 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); + boundary_condition tmp_bc = bc(end + 1, j + 1); if (tmp_bc.type == Diffusion::BC_CLOSED) continue; @@ -224,7 +224,8 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c, y_values[2] = 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); + (y_values[0] - 3 * y_values[1] + 2 * y_values[2]) / + (dx * dx); } return t0_c; From 437ab1b10b65ee80196c012d1f9562ab6f8ccd24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 27 Apr 2022 12:50:37 +0200 Subject: [PATCH 34/36] Update .gitlab-ci.yml file --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 2e36801..872330f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -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 From 56e70efd1ebe092b97656515c3f4c014d3e90c71 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Wed, 27 Apr 2022 12:54:01 +0200 Subject: [PATCH 35/36] Update CI file to build apps with debug symbols --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 872330f..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: From b2bbd281756111ce88c47b545a429d8847ea90d4 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 27 Apr 2022 14:12:43 +0200 Subject: [PATCH 36/36] Fix clang-tidy suggestions --- include/diffusion/BTCSDiffusion.hpp | 19 ++++++++++--------- src/BTCSDiffusion.cpp | 8 +++++--- 2 files changed, 15 insertions(+), 12 deletions(-) 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 54b9789..240d76a 100644 --- a/src/BTCSDiffusion.cpp +++ b/src/BTCSDiffusion.cpp @@ -179,14 +179,15 @@ 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++) { boundary_condition tmp_bc = bc(0, j + 1); - if (tmp_bc.type == Diffusion::BC_CLOSED) + 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); @@ -216,8 +217,9 @@ auto Diffusion::BTCSDiffusion::calc_t0_c(const DMatrixRowMajor &c, for (int j = 0; j < n_cols; j++) { boundary_condition tmp_bc = bc(end + 1, j + 1); - if (tmp_bc.type == Diffusion::BC_CLOSED) + if (tmp_bc.type == Diffusion::BC_CLOSED) { continue; + } y_values[0] = c(end - 1, j); y_values[1] = c(end, j);