Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix rendering docs to PDF #737

Merged
merged 4 commits into from
Nov 6, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/core/allocation.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -172,7 +172,7 @@ $$
$$
\hat{C}_S^p \leftarrow \hat{C}_S^p - F;
$$
5. Repeat steps 2-4 for the remaining priorities up to $p_\max$.
5. Repeat steps 2-4 for the remaining priorities up to $p_{\max}$.

:::{.callout-note}
In the future there will be another solve before the priority 1 solve, taking the demand/supply from basins into account.
Expand Down
34 changes: 15 additions & 19 deletions docs/core/equations.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -94,18 +94,16 @@ The presence of division by the basin area means that areas of size zero are not
## The reduction factor {#sec-reduction_factor}
At several points in the equations below a *reduction factor* is used. This is a term that makes certain transitions more smooth, for instance when a pump stops providing water when its source basin dries up. The reduction factor is given by

$$
\begin{align}
\phi(x; p) =
\begin{align}
\begin{cases}
0 &\text{if}\quad x < 0 \\
-2 \left(\frac{x}{p}\right)^3 + 3\left(\frac{x}{p}\right)^2 &\text{if}\quad 0 \le x \le p \\
1 &\text{if}\quad x > p
\end{cases}
\end{align},
$$
\begin{cases}
0 &\text{if}\quad x < 0 \\
-2 \left(\frac{x}{p}\right)^3 + 3\left(\frac{x}{p}\right)^2 &\text{if}\quad 0 \le x \le p \\
1 &\text{if}\quad x > p
\end{cases}
\end{align}

where $p > 0$ is the threshold value which determines the interval $[0,p]$ of the smooth transition between $0$ and $1$, see the plot below.
Here $p > 0$ is the threshold value which determines the interval $[0,p]$ of the smooth transition between $0$ and $1$, see the plot below.

```{python}
# | code-fold: true
Expand Down Expand Up @@ -463,23 +461,21 @@ $$
Q_\text{PID}(t) = K_p e(t) + K_i\int_{t_0}^t e(\tau)\text{d}\tau + K_d \frac{\text{d}e}{\text{d}t},
$$ {#eq-PIDflow}

for given constant parameters $K_p,K_i,K_d$. The pump or outlet can have associated minimum and maximum flow rates $Q_\min, Q_\max$, and so
for given constant parameters $K_p,K_i,K_d$. The pump or outlet can have associated minimum and maximum flow rates $Q_{\min}, Q_{\max}$, and so
$$
Q_\text{pump/outlet} = \text{clip}(\Phi Q_\text{PID}; Q_\min, Q_\max).
Q_\text{pump/outlet} = \text{clip}(\Phi Q_\text{PID}; Q_{\min}, Q_{\max}).
$$

Here $u_\text{us}$ is the storage of the basin upstream of the pump or outlet, $\Phi$ is the product of [reduction factors](equations.qmd#sec-reduction_factor) associated with the [pump](equations.qmd#sec-pump) or [outlet](equations.qmd#sec-outlet) and

$$
\text{clip}(Q; Q_\min, Q_\max) =
\begin{align}
\text{clip}(Q; Q_{\min}, Q_{\max}) =
\begin{cases}
Q_\min & \text{if} & Q < Q_\min \\
Q & \text{if} & Q_\min \leq Q \leq Q_\max \\
Q_\max & \text{if} & Q > Q_\max
Q_{\min} & \text{if} \quad Q < Q_{\min} \\
Q & \text{if} \quad Q_{\min} \leq Q \leq Q_{\max} \\
Q_{\max} & \text{if} \quad Q > Q_{\max}
\end{cases}.
\end{align}
$$

For the integral term we denote
$$
Expand Down Expand Up @@ -511,7 +507,7 @@ that is, $\hat{f}_\text{PID}$ is the right hand side of the ODE for the controll

Using this, solving @eq-PIDflow for $Q_\text{PID}$ yields
$$
Q_\text{pump/outlet} = \text{clip}\left(\phi(u_\text{us})\frac{K_pe + K_iI + K_d \left(\frac{\text{d}\text{SP}}{\text{d}t}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right)}{1\pm\phi(u_\text{us})\frac{K_d}{A(u_\text{PID})}};Q_\min,Q_\max\right),
Q_\text{pump/outlet} = \text{clip}\left(\phi(u_\text{us})\frac{K_pe + K_iI + K_d \left(\frac{\text{d}\text{SP}}{\text{d}t}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right)}{1\pm\phi(u_\text{us})\frac{K_d}{A(u_\text{PID})}};Q_{\min},Q_{\max}\right),
$$
where the clipping is again done last. Note that to compute this, $\hat{f}_\text{PID}$ has to be known first, meaning that the PID controlled pump/outlet flow rate has to be computed after all other contributions to the PID controlled basin's storage are known.

Expand Down
9 changes: 4 additions & 5 deletions docs/core/numerics.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -49,25 +49,24 @@ $$
\mathbf{g}(\mathbf{x}) = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{x},t_{n+1})
$$
by iterating $\mathbf{g}$ on an initial guess of $\mathbf{w}_{n+1}$;
- Newton-Raphson iterations: approximate $\mathbf{w}_{n+1}$ as a root of the function
- Newton iterations: approximate $\mathbf{w}_{n+1}$ as a root of the function
$$
\mathbf{h}(\mathbf{x}) = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{x},t_{n+1}) - \mathbf{x},
$$
by iteratively finding the root of its linearized form:
$$

\begin{align}
\mathbf{0} =& \mathbf{h}(\mathbf{w}_{n+1}^k) + \mathbf{J}(\mathbf{h})(\mathbf{w}_{n+1}^k)(\mathbf{w}_{n+1}^{k+1}-\mathbf{w}_{n+1}^k) \\
=& \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{w}_{n+1}^k,t_{n+1}) - \mathbf{w}_{n+1}^k \\ +&\left[(t_{n+1}-t_n)\mathbf{J}(\mathbf{f})(\mathbf{w}_{n+1}^k)-\mathbf{I}\right](\mathbf{w}_{n+1}^{k+1}-\mathbf{w}_{n+1}^k).
\end{align}
$$ {#eq-newtoniter}
Note that this thus requires an evaluation of the Jacobian of $\mathbf{f}$ and solving a linear system per iteration.

# The advantage of implicit methods

The implicit method @eq-eulerbackward is a coupled system of equations for $\mathbf{w}_{n+1}$, while the explicit @eq-eulerforward is fully decoupled. This means in general that in a time integration step with an implicit method the basins communicate information with eachother, while in an explicit method this does not happen. A consequence of this is that local events (e.g. a pump turns on) propagate slowly trough the model using an explicit method but quickly using an implicit method, making implicit methods more stable.

# Jacobian computations
The iterations @eq-newtoniter require an evaluation of the Jacobian of $\mathbf{f}$. The Jacobian of `water_balance!` is discussed [here](equations.qmd#the-jacobian).
The Newton iterations above require an evaluation of the Jacobian of $\mathbf{f}$. The Jacobian of `water_balance!` is discussed [here](equations.qmd#the-jacobian).

There are several ways to compute the Jacobian:

Expand All @@ -77,7 +76,7 @@ There are several ways to compute the Jacobian:

# Continuity considerations

The convergence of the Newton-Raphson method can be [proven](https://en.wikipedia.org/wiki/Newton%27s_method#Proof_of_quadratic_convergence_for_Newton's_iterative_method) given certain properties of $\mathbf{f}$ around the initial guess and root to find. An important aspect is the smoothness of $\mathbf{f}$. The basin profiles and $Q(h)$ relations are given by interpolated data, and thus we have some control over the smoothness of these functions by the choice of interpolation method. This is discussed further below. the Manning resistance is not mentioned here since it is given by an analytical expression.
The convergence of the Newton method can be [proven](https://en.wikipedia.org/wiki/Newton%27s_method#Proof_of_quadratic_convergence_for_Newton's_iterative_method) given certain properties of $\mathbf{f}$ around the initial guess and root to find. An important aspect is the smoothness of $\mathbf{f}$. The basin profiles and $Q(h)$ relations are given by interpolated data, and thus we have some control over the smoothness of these functions by the choice of interpolation method. This is discussed further below. the Manning resistance is not mentioned here since it is given by an analytical expression.

Control mechanisms can change parameters of $\mathbf{f}$ discontinuously, leading to discontinuities of $\mathbf{f}$. This however does not yield problems for the time integration methods in `DifferentialEquations.jl`, since the [callback mechanisms](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) used to change these parameters make the solver take these discontinuities into account.

Expand Down
8 changes: 0 additions & 8 deletions docs/couple/modflow.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -235,12 +235,10 @@ groundwater response to Ribasim.
From Ribasim's perspective, the groundwater head is constant given a timestep,
so that:

$$
\begin{align}
p = -C \\
q = -C\phi
\end{align}
$$

When the head falls below the drainage elevation, the coefficients are 0.

Expand All @@ -249,22 +247,18 @@ When the head falls below the drainage elevation, the coefficients are 0.
From Ribasim's perspective, infiltration is never limited when the head falls below
the bottom:

$$
\begin{align}
p = -C \\
q = -Cb
\end{align}
$$

Otherwise, infiltration and drainage occur with the same equation as for the
drainage package:

$$
\begin{align}
p = -C \\
q = -C\phi
\end{align}
$$


# Assigning MODFLOW water levels
Expand All @@ -273,12 +267,10 @@ Principally, MODFLOW deals with water levels; Ribasim deals with water volumes.
The translation from a water volume to a water level for a single boundary
condition is straightforward:

$$
\begin{align}
\Delta d = \frac{\Delta V}{A} \\
h = b + d
\end{align}
$$

Where $d$ is water depth, $V$ is volume, $A$ is area, $h$ is water level, and
$b$ is the bed elevation.
Expand Down