-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'master' of https://github.com/storm-fsv-cvut/smoderp2d-…
- Loading branch information
Showing
7 changed files
with
470 additions
and
468 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,3 +1,236 @@ | ||
| \subsection{Source code architecture} | ||
|
|
||
| \subsection{NumPy solution} | ||
| \subsection{Implicit method numerical solution} | ||
|
|
||
| \subsection{Explicit/Implicit approach} | ||
|
|
||
| \subsubsection{Explicit} | ||
| The time derivative in Equation \ref{eq:bilance} is calculated using an | ||
| explicit method. The computation is therefore sensitive to the size of the time | ||
| step. The size of the time step is controlled by the Courant criterion, which | ||
| needs to be kept below the theoretical maximum value of 1, while the maximum | ||
| value in practise is lower than 1 | ||
| \cite{zhang1989modeling, esteves2000overland}. | ||
|
|
||
|
|
||
| The sheet flow water level of the next time t + 1 step in Equation | ||
| \ref{eq:bilance} which incorporates the sum \ref{eq:routing} is calculated with the | ||
| explicit time discretisation scheme for cell i as: | ||
| \begin{equation} | ||
| h_{i,t} =h_{i,t-1} + \mathrm{d}t (es_{i,t-1} + \sum_j^m q^{out}_{j,t-1}- | ||
| inf_{i,t-1} - q^{out}_{i,t-1}), | ||
| \label{eq:bilanceexpl} | ||
| \end{equation} | ||
|
|
||
|
|
||
| After inserting the equations~(\ref{eq:infiltration}) | ||
| and~(\ref{eq:powerlaw2}) into equation~(\ref{eq:bilanceexpl}) the | ||
| final explicit form balance equation can be written as: | ||
| \begin{dmath} | ||
| h_{i,t} = h_{i,t-1} + | ||
| es_{i,t-1} + exf_{i,t-1} + \sum_{j}^{inflows} | ||
| 1/n_{sheet,j}I_j^{y_j} h_{j,t-1}^{b_j} - | ||
| 1/n_{sheet,i}I_i^{y_i}h_{i,t-1}^{b_i} - ( | ||
| \frac{1}{2}S_it_1^{-1/2}+K_{s,i}). | ||
| \end{dmath} | ||
|
|
||
|
|
||
| \subsubsection{Implicit} | ||
| Alternatively the equation (\ref{eq:bilance}) solved using implicit | ||
| scheme. The right hand side can be expressed in for time $t$ as: | ||
|
|
||
| $$ | ||
| \frac{h_{i,t} - h_{i,t-1} }{\mathrm{d} t} = \left(es_{i,t} + | ||
| \sum_j^{inflows} q^{in}_{j,t} - inf_{i,t} - q^{out}_{i,t}\right). | ||
| $$ | ||
|
|
||
| After inserting the power-law equation (\ref{eq:powerlaw}) and several rearrangements the formula above can be written | ||
| as follows: | ||
|
|
||
|
|
||
|
|
||
| $$ | ||
| h_{i,t} + \mathrm{d} t a_ih^{b_{i}-1}_{i,t} h_{i,t} - \mathrm{d} t \sum_j^{inflows} | ||
| a_jh^{b_{j}-1}_{j,t} h_{j,t} = h_{i,t-1} + \mathrm{d} t es_{i,t} - \mathrm{d} | ||
| t inf_{i,t}. | ||
| $$ | ||
| Here the known part of the balance equation (rainfall and infiltration) are at the right hand side. To | ||
| extract the unknown $h$ form the power law a following rearrangement was used: | ||
|
|
||
|
|
||
|
|
||
| \begin{equation} | ||
| (1+\mathrm{d} t a_ih^{b_{i}-1}_{i,t})h_{i,t} - \mathrm{d} t \sum_j^{inflows} | ||
| a_jh^{b_{j}-1}_{j,t} h_{j,t} = h_{i,t-1} + \mathrm{d} t es_{i,t} - | ||
| \mathrm{d} t inf_{i,t} | ||
| \label{eq:impl} | ||
| \end{equation} | ||
|
|
||
| From the equation (\ref{eq:impl}) the system of linear equations can be constructed. | ||
|
|
||
| As example, lets show the equation (\ref{eq:impl}) for | ||
| $ | ||
| i=5\ \mathrm{and}\ inflows\in\{7,8,9\} | ||
| $: | ||
|
|
||
|
|
||
| \begin{dmath} | ||
| (1+\Delta T a_5h^{b_{5}-1}_{5,t})h_{5,t} - \Delta T a_7h^{b_{7}-1}_{7,t} h_{7,t} + \Delta T a_8h^{b_{8}-1}_{8,t} h_{8,t} + \Delta T a_9 h^{b_{9}-1}_{9,t} h_{9,t} = h_{5,t-1} + \Delta T es_{5,t} - \Delta T inf_{5,t}. | ||
| \end{dmath} | ||
|
|
||
| In matrix form the above equation can be written as: | ||
|
|
||
|
|
||
|
|
||
| \begin{dmath} | ||
| \begin{bmatrix} | ||
| \hdots& & & & & & \\ | ||
| \hdots & 1+\Delta T a_5h^{b_{5}-1}_{5,t} & \hdots & \Delta T a_7h^{b_{7}-1}_{7,t} & \Delta T a_8h^{b_{8}-1}_{8,t} & \Delta T a_9 h^{b_{9}-1}_{9,t} & \hdots \\ | ||
| \hdots& & & & & & \\ | ||
| \hdots& && & & & \\ | ||
| \hdots& && & & & \\ | ||
| \hdots& && & & & \\ | ||
|
|
||
| \end{bmatrix} | ||
| % | ||
| \begin{bmatrix} | ||
| \vdots \\ | ||
| h_{5,t} \\ | ||
| \vdots \\ | ||
| h_{7,t} \\ | ||
| h_{8,t} \\ | ||
| h_{9,t} \\ | ||
| \vdots \\ | ||
| % | ||
| \end{bmatrix} | ||
| = | ||
| \begin{bmatrix} | ||
| \vdots \\ | ||
| h_{5,t-1} + \Delta T es_{5,t} - \Delta T inf_{5,t} \\ | ||
| \vdots \\ | ||
| \vdots \\ | ||
| \vdots \\ | ||
| \vdots \\ | ||
| \vdots \\ | ||
| % | ||
| \end{bmatrix} | ||
| \label{eq:sys} | ||
| \end{dmath} | ||
|
|
||
| The system (\ref{eq:sys}) is solved using \texttt{scipy} package method \texttt{root} | ||
| using \texttt{krigov} methof that shown the best convergence. | ||
|
|
||
|
|
||
| \paragraph{Implicit solution with rill flow} To calculated rill flow the water level in rill and needs to be calculated as | ||
| \begin{equation} | ||
| h_{rill} = \text{max}(h-h_{crit},0), | ||
| \label{eq:hrill2} | ||
| \end{equation} | ||
| where | ||
| \begin{equation} | ||
| h_{sheet} = \text{min}(h,h_{crit}). | ||
| \label{eq:hsheet2} | ||
| \end{equation} | ||
|
|
||
| Mannings formula is used to calculated the steady state flow in | ||
| the rill assuming rill is a channel (description in section \ref{sec:rill}. | ||
| $$ | ||
| q_{rill} = 1/n R(h_{rill})^{2/3} i^{1/2} | ||
| $$ | ||
|
|
||
|
|
||
| The balance equation becomes | ||
| \begin{dmath} | ||
| \frac{h_{i,t} - h_{i,t-1} }{\Delta T} = | ||
| es_{i,t} + \sum_j^{inflows} a_jh^{b_{j}}_{sh,j,t} + \sum_k^{inflows} 1/n_k R_k(h_{rl,k,t})^{2/3} i_k^{1/2} - inf_{i,t} - a_ih^{b_{i}}_{sh,i,t} - 1/n R(h_{rl,i,t})^{2/3} i^{1/2}, | ||
| \end{dmath} | ||
| where $h_{sheet}$ and $h_{rill}$ are defined as (\ref{eq:hsheet2}) and (\ref{eq:hrill2}). | ||
|
|
||
| This formula is rearranged and, for practical purposes, the linear equation is the system is calculated with condition \\ | ||
| for $h<=h_{crit}$ as | ||
| \begin{equation} | ||
| \left(\frac{1}{\Delta T}+a_ih^{b_{i}-1}_{i,t}\right)h_{i,t} - \sum_j^{inflows} a_jh^{b_{j}-1}_{j,t} h_{j,t} = \frac{h_{i,t-1}}{\Delta} + es_{i,t} - inf_{i,t} | ||
| \end{equation} | ||
| % | ||
| % | ||
| and for $h>h_{crit}$ as | ||
| \begin{multline} | ||
| \left(\frac{1}{\Delta T} | ||
| + 1/n_i R_i(h_{i,t}-h_{crit,i})^{2/3} i_i^{1/2} \frac{1}{h_{i,t}}\right)h_{i,t} \\ | ||
| - \sum_k^{inflows} \left( 1/n_k R_k(h_{k,t}-h_{crit,k})^{2/3} i_k^{1/2} \frac{1}{h_{k,t}}\right)h_{k,t} | ||
| = \\ | ||
| \frac{h_{i,t-1} }{\Delta T} | ||
| + es_{i,t} | ||
| + \sum_j^{inflows} a_j h^{b_{j}}_{crit,j} | ||
| - inf_{i,t} | ||
| - a_i h^{b_{i}}_{crit,i} | ||
| \end{multline} | ||
| The later system of equations is constructed in similar fashion as the case of (\ref{eq:sys}) and solved with the same solver. | ||
|
|
||
|
|
||
|
|
||
| \subsection{Solving instability in calculation} | ||
|
|
||
| At higher flow velocities or with inappropriately chosen time step or pixel | ||
| size, or if abrupt change in solution occurs, instabilities occurred in | ||
| the solution. In the current version of the SMODERP2D program, this problem | ||
| is resolved for explicit solution by the Courant-Friedrich-Lewy (CFL) | ||
| condition and for implicit solution the model change the time step based on | ||
| the number of iterations in the system solution. | ||
|
|
||
|
|
||
| \subsubsection{Courant-Friedrich-Lewy condition} | ||
| Filling of this condition ensures convergence of the explicit solution if | ||
| $\text{CFL} < 1.0$. The general equation of the CFL condition was derived | ||
| and adapted for the purposes of the SMODERP2D model to the following form: | ||
|
|
||
| \begin{equation} | ||
| \text{CFL} = \frac{1}{0.5601}\frac{v \Delta T}{\Delta X} | ||
| \label{eq:courrovnice} | ||
| \end{equation} where $v$ represents the velocity of surface or channel flow | ||
| [$m/s$], $\Delta T$ and $\Delta X$ represent time step and spatial step | ||
| respectively. The number 0.5601 was empiricaly shown to improve the | ||
| stability. Lower than 1 CFL conditions was also indicated by | ||
| \cite{zhang1989modeling, esteves2000overland}. | ||
|
|
||
| After completing of single time step, the highest value of CFL determined | ||
| from surface runoff using equation (\ref{eq:courrovnice}) is stored. Then, | ||
| this value is compared to the critical value of CFL, and according to the | ||
| rules illustrated in Table \ref{tab:cflsheet}, the time step length $\Delta | ||
| T$ is changed (or not changed). If $\Delta T$ changes, the calculation is | ||
| repeated in the given time step using the updated $\Delta T$. The | ||
| calculation progresses to the next time only when the stability of the | ||
| calculation is guaranteed. | ||
|
|
||
| \begin{table}[t!] | ||
| \centering | ||
| \caption{Criteria for changing the time step based on surface runoff} | ||
| \label{tab:cflsheet} | ||
| \begin{tabular}{ccc} | ||
| \hline | ||
| New & $\text{CFL} < 0.75 \lor 1.0 < \text{CFL}$ & $ 0.75 \geq \text{CFL} \geq 1.0 \lor \text{CFL} = 0.0^*$ \\ | ||
| \hline | ||
| \hline | ||
| $\Delta T$ & = $\min\left(\frac{0.5601 \Delta X}{v}, \Delta T_{\text{max}}\right)$ & = original $\Delta T$\\ | ||
| \hline | ||
| % \multicolumn{3}{l}{{\small CFL = 0.0 typically in the case when the flow velocity is zero. In that case,}} | ||
| \end{tabular} | ||
| \end{table} | ||
|
|
||
|
|
||
| \subsubsection{Stability in implicit solution} | ||
| To ensure the stability in implicit solution, the model checks the number | ||
| of iterations in the in the solver of the linear system. High number of | ||
| iteration indicate abrupt change in the solution, that leads to instability | ||
| in the solution. The changes in time step length follow simple rules, that | ||
| are hard-coded to the model in the current version. The rules are as | ||
| followsL | ||
| \begin{itemize} | ||
| \item if the number of iterations exceeds $it_{max}$, then the time step | ||
| is multiplied by factor smaller then one. | ||
| \item if the number of iterations is lower then $it_{min}$, then the time | ||
| step is multiplied by factor larger then one. | ||
| \item if the number of iterations is higher then $it_{min}$ and lower | ||
| then $it_{max}$, then the time step remains unchanged. | ||
| \end{itemize} | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.