deploy: 221d51b

concept/allocation.html

@@ -578,7 +578,7 @@ <h2 data-number="4.3" class="anchored" data-anchor-id="the-optimization-constrai
 <section id="example" class="level2" data-number="4.4">
 <h2 data-number="4.4" class="anchored" data-anchor-id="example"><span class="header-section-number">4.4</span> Example</h2>
 <p>The following is an example of an optimization problem for the example shown <a href="../guide/examples.html#model-with-allocation-user-demand">here</a>:</p>
-<div id="8dd273bf" class="cell" data-execution_count="1">
+<div id="b34ec3e8" class="cell" data-execution_count="1">
 <details class="code-fold">
 <summary>Code</summary>
 <div class="sourceCode cell-code" id="cb1"><pre class="sourceCode julia code-with-copy"><code class="sourceCode julia"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="im">using</span> <span class="bu">Ribasim</span></span>
@@ -602,25 +602,25 @@ <h2 data-number="4.4" class="anchored" data-anchor-id="example"><span class="hea
 <div class="cell-output cell-output-stdout">
 <pre><code>Min F[(Basin #2, UserDemand #3)]² + F[(Basin #5, UserDemand #6)]²
 Subject to
- flow_conservation[Basin #5] : -F[(Basin #5, TabulatedRatingCurve #7)] + F[(LinearResistance #4, Basin #5)] - F[(Basin #5, LinearResistance #4)] + F[(UserDemand #6, Basin #5)] - F[(Basin #5, UserDemand #6)] = 0
- flow_conservation[Terminal #8] : F[(TabulatedRatingCurve #7, Terminal #8)] = 0
- flow_conservation[LinearResistance #4] : -F[(LinearResistance #4, Basin #5)] + F[(Basin #5, LinearResistance #4)] + F[(Basin #2, LinearResistance #4)] - F[(LinearResistance #4, Basin #2)] = 0
- flow_conservation[Basin #2] : F[(FlowBoundary #1, Basin #2)] - F[(Basin #2, UserDemand #3)] + F[(UserDemand #3, Basin #2)] - F[(Basin #2, LinearResistance #4)] + F[(LinearResistance #4, Basin #2)] = 0
  flow_conservation[TabulatedRatingCurve #7] : F[(Basin #5, TabulatedRatingCurve #7)] - F[(TabulatedRatingCurve #7, Terminal #8)] = 0
+ flow_conservation[Basin #2] : -F[(Basin #2, UserDemand #3)] + F[(UserDemand #3, Basin #2)] + F[(FlowBoundary #1, Basin #2)] - F[(Basin #2, LinearResistance #4)] + F[(LinearResistance #4, Basin #2)] = 0
+ flow_conservation[Basin #5] : -F[(Basin #5, TabulatedRatingCurve #7)] + F[(LinearResistance #4, Basin #5)] - F[(Basin #5, LinearResistance #4)] - F[(Basin #5, UserDemand #6)] + F[(UserDemand #6, Basin #5)] = 0
+ flow_conservation[LinearResistance #4] : -F[(LinearResistance #4, Basin #5)] + F[(Basin #5, LinearResistance #4)] + F[(Basin #2, LinearResistance #4)] - F[(LinearResistance #4, Basin #2)] = 0
+ flow_conservation[Terminal #8] : F[(TabulatedRatingCurve #7, Terminal #8)] = 0
  source[(FlowBoundary #1, Basin #2)] : F[(FlowBoundary #1, Basin #2)] ≤ 172800
  source_user[UserDemand #6] : F[(UserDemand #6, Basin #5)] ≤ 0
  source_user[UserDemand #3] : F[(UserDemand #3, Basin #2)] ≤ 0
- F[(FlowBoundary #1, Basin #2)] ≥ 0
+ F[(Basin #2, UserDemand #3)] ≥ 0
  F[(Basin #5, TabulatedRatingCurve #7)] ≥ 0
+ F[(UserDemand #3, Basin #2)] ≥ 0
  F[(LinearResistance #4, Basin #5)] ≥ 0
  F[(Basin #5, LinearResistance #4)] ≥ 0
- F[(UserDemand #6, Basin #5)] ≥ 0
- F[(Basin #2, UserDemand #3)] ≥ 0
- F[(TabulatedRatingCurve #7, Terminal #8)] ≥ 0
  F[(Basin #5, UserDemand #6)] ≥ 0
- F[(UserDemand #3, Basin #2)] ≥ 0
+ F[(FlowBoundary #1, Basin #2)] ≥ 0
+ F[(TabulatedRatingCurve #7, Terminal #8)] ≥ 0
  F[(Basin #2, LinearResistance #4)] ≥ 0
  F[(LinearResistance #4, Basin #2)] ≥ 0
+ F[(UserDemand #6, Basin #5)] ≥ 0
 </code></pre>
 </div>
 </div>

concept/equations.html

@@ -269,8 +269,12 @@ <h2 id="toc-title">On this page</h2>
   <li><a href="#the-pid-control-integral-state" id="toc-the-pid-control-integral-state" class="nav-link" data-scroll-target="#the-pid-control-integral-state"><span class="header-section-number">1.1</span> The PID control integral state</a></li>
   <li><a href="#exactly-integrating-flows-to-minimize-the-number-of-states" id="toc-exactly-integrating-flows-to-minimize-the-number-of-states" class="nav-link" data-scroll-target="#exactly-integrating-flows-to-minimize-the-number-of-states"><span class="header-section-number">1.2</span> Exactly integrating flows to minimize the number of states</a></li>
   <li><a href="#the-jacobian" id="toc-the-jacobian" class="nav-link" data-scroll-target="#the-jacobian"><span class="header-section-number">1.3</span> The Jacobian</a></li>
-  <li><a href="#why-this-formulation" id="toc-why-this-formulation" class="nav-link" data-scroll-target="#why-this-formulation"><span class="header-section-number">1.4</span> Why this formulation</a></li>
-  <li><a href="#numerical-solution" id="toc-numerical-solution" class="nav-link" data-scroll-target="#numerical-solution"><span class="header-section-number">1.5</span> Numerical solution</a></li>
+  <li><a href="#the-water-balance-error" id="toc-the-water-balance-error" class="nav-link" data-scroll-target="#the-water-balance-error"><span class="header-section-number">1.4</span> The water balance error</a>
+  <ul class="collapse">
+  <li><a href="#example-calculation" id="toc-example-calculation" class="nav-link" data-scroll-target="#example-calculation"><span class="header-section-number">1.4.1</span> Example calculation</a></li>
+  </ul></li>
+  <li><a href="#why-this-formulation" id="toc-why-this-formulation" class="nav-link" data-scroll-target="#why-this-formulation"><span class="header-section-number">1.5</span> Why this formulation</a></li>
+  <li><a href="#numerical-solution" id="toc-numerical-solution" class="nav-link" data-scroll-target="#numerical-solution"><span class="header-section-number">1.6</span> Numerical solution</a></li>
   </ul></li>
   <li><a href="#performance" id="toc-performance" class="nav-link" data-scroll-target="#performance"><span class="header-section-number">2</span> Performance</a></li>
   </ul>
@@ -346,13 +350,63 @@ <h2 data-number="1.3" class="anchored" data-anchor-id="the-jacobian"><span class
     J_{i,j} = \frac{\partial f_j}{\partial u_i},
 \]</span> i.e.&nbsp;<span class="math inline">\(J_{i,j}\)</span> quantifies how <span class="math inline">\(f_j\)</span>. the time derivative of state <span class="math inline">\(j\)</span>, changes with respect to changes in state <span class="math inline">\(i\)</span>. Most of these entries are <span class="math inline">\(0\)</span>, because flows in distant parts of the model do not depend on each other.</p>
 </section>
-<section id="why-this-formulation" class="level2" data-number="1.4">
-<h2 data-number="1.4" class="anchored" data-anchor-id="why-this-formulation"><span class="header-section-number">1.4</span> Why this formulation</h2>
+<section id="the-water-balance-error" class="level2" data-number="1.4">
+<h2 data-number="1.4" class="anchored" data-anchor-id="the-water-balance-error"><span class="header-section-number">1.4</span> The water balance error</h2>
+<p>The water balance error quantifies how well the water volume in the model is conserved for each Basin over an output save period, i.e.&nbsp;whether no water erroneously appears or disappears. It looks at the storage rate <span class="math display">\[
+    \text{storage rate} = \frac{\Delta S_b}{\Delta t}
+\]</span></p>
+<p>in a Basin over a time period <span class="math inline">\(\Delta t = t_1 - t_0\)</span> and compares that to the total inflows and outflows of that Basin over that period. More precisely, we first compute the total inflow and outflow, where:</p>
+<ul>
+<li><span class="math inline">\(\text{total inflow}\)</span>: the precipitation, drainage and horizontal flows into the Basin;</li>
+<li><span class="math inline">\(\text{total outflow}\)</span>: the evaporation, infiltration and horizontal flows out of the Basin.</li>
+</ul>
+<p>Whether a flow is an inflow or an outflow depends on whether the flow contributes to or takes from the Basin storage, and thus this is independent of the edge direction. This is determined for each solver timestep individually.</p>
+<p>Then from this we compute the errors:</p>
+<p><span class="math display">\[
+    \begin{align}
+    \text{balance error} =&amp;&amp; \text{storage rate} - (\text{total inflow} - \text{total outflow}) \\
+    \text{relative error}=&amp;&amp; \frac{\text{absolute error}}{0.5(\text{total inflow} + \text{total outflow})}
+    \end{align}
+\]</span> Hence the reference used for computing the relative error is the average of the total inflow and total outflow of the Basin (which are both non-negative).</p>
+<p>The default tolerances are <span class="math inline">\(0.001 \text{ m}^3\)</span> for the balance error and <span class="math inline">\(0.01\)</span> for the absolute error, which should not be exceeded for realistic models.</p>
+<p>In extreme cases where storage rate is many orders of magnitude smaller than the storage itself, these computations can have floating point truncation errors which can lead to large relative errors. This is however only when the storage is roughly <span class="math inline">\(10^{15}\)</span> or more times bigger than the storage rate.</p>
+<section id="example-calculation" class="level3" data-number="1.4.1">
+<h3 data-number="1.4.1" class="anchored" data-anchor-id="example-calculation"><span class="header-section-number">1.4.1</span> Example calculation</h3>
+<p>Say we have the following model:</p>
+<div class="quarto-figure quarto-figure-left">
+<figure class="figure">
+<p><img src="https://s3.deltares.nl/ribasim/doc-image/concept/equations/subnetwork.png" class="img-fluid quarto-figure quarto-figure-left figure-img"></p>
+</figure>
+</div>
+<p>and we want to calculate the water balance error for Basin 6. We have the following data:</p>
+<ul>
+<li>Time period length: <span class="math inline">\(10.0 \text{ s}\)</span></li>
+<li>Basin storage start: <span class="math inline">\(100.0 \text{ m}^3\)</span></li>
+<li>Basin storage end: <span class="math inline">\(50.0 \text{ m}^3\)</span></li>
+<li>UserDemand #11 inflow average: <span class="math inline">\(10.0 \text{ m}^3/\text{s}\)</span></li>
+<li>UserDemand #11 outflow average: <span class="math inline">\(5.0 \text{ m}^3/\text{s}\)</span></li>
+<li>Outlet #7 flow average: <span class="math inline">\(- 3.5 \text{ m}^3/\text{s}\)</span></li>
+<li>Outlet #11 flow average: <span class="math inline">\(4.0 \text{ m}^3/\text{s}\)</span></li>
+</ul>
+<p>And so we get</p>
+<p><span class="math display">\[
+\begin{align}
+    \text{storage rate} = &amp;&amp; \frac{50.0 - 100.0}{10.0} &amp;= &amp; -6.0 \text{ m}^3/\text{s} \\
+    \text{total inflow} = &amp;&amp; 5.0 + 3.5 &amp;= &amp; 8.5 \text{ m}^3/\text{s}\\
+    \text{total outflow} = &amp;&amp; 10.0 + 4.0 &amp;= &amp; 14.0 \text{ m}^3/\text{s}\\
+    \text{balance error} = &amp;&amp; -6.0 - (8.5 - 14.0) &amp;= &amp; -0.5 \text{ m}^3/\text{s}\\
+    \text{relative error} = &amp;&amp; \frac{-0.5}{8.5 + 14.0} &amp;\approx &amp; -0.022
+\end{align}
+\]</span> Note that the balance error and relative error are negative, but we use their absolute value to compare to the respective tolerances.</p>
+</section>
+</section>
+<section id="why-this-formulation" class="level2" data-number="1.5">
+<h2 data-number="1.5" class="anchored" data-anchor-id="why-this-formulation"><span class="header-section-number">1.5</span> Why this formulation</h2>
 <p>You might wonder why in the above explanation the states are given by the cumulative flows and not by the Basin storages, which is arguably conceptually simpler. The reason is that we do not just want to model the storages in the Basins over time, but we also want accurate output of each individual flow, e.g.&nbsp;to model the spread of pollutants.</p>
 <p>When the states are given by the storages, generally the individual flows can not accurately be computed from that as a post processing step, because there are more flows than storages. Also, we can only compute flows at individual points in time explicitly, not over a whole interval. When the states are given by the cumulative flows however, the output of the problem solve gives these flows directly, and from those the storage over time can be computed accurately. Hence in short, the formulation above gives more information than a formulation with Basin storages as states.</p>
 </section>
-<section id="numerical-solution" class="level2" data-number="1.5">
-<h2 data-number="1.5" class="anchored" data-anchor-id="numerical-solution"><span class="header-section-number">1.5</span> Numerical solution</h2>
+<section id="numerical-solution" class="level2" data-number="1.6">
+<h2 data-number="1.6" class="anchored" data-anchor-id="numerical-solution"><span class="header-section-number">1.6</span> Numerical solution</h2>
 <p>Ribasim uses <a href="https://github.com/SciML/OrdinaryDiffEq.jl/">OrdinaryDiffEq.jl</a> to provide a numerical solution to the water balance equations. Changes to forcings or parameters such as precipitation, but also the allocated water abstraction is managed through the use of callback functions <span class="citation" data-cites="callbacks">(<a href="#ref-callbacks" role="doc-biblioref">SciML Development Team 2022</a>)</span>. In a coupled run, the exchanges with MODFLOW 6 are also managed via the use of a callback function. For more a more in-depth discussion of numerical computations see <a href="../concept/numerics.html">Numerical considerations</a>.</p>
 </section>
 </section>

reference/node/discrete-control.html

@@ -383,7 +383,7 @@ <h1 class="title">DiscreteControl</h1>
 
 
 <p>Set parameters of other nodes based on model state conditions (e.g.&nbsp;Basin level). The table below shows which parameters are controllable for a given node type.</p>
-<div id="00be1382" class="cell" data-execution_count="1">
+<div id="6acb3fb6" class="cell" data-execution_count="1">
 <details class="code-fold">
 <summary>Code</summary>
 <div class="sourceCode cell-code" id="cb1"><pre class="sourceCode julia code-with-copy"><code class="sourceCode julia"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="im">using</span> <span class="bu">Ribasim</span></span>

reference/node/manning-resistance.html

@@ -539,7 +539,7 @@ <h1 data-number="2"><span class="header-section-number">2</span> Equations</h1>
   \end{cases}
 \end{align}
 \]</span></p>
-<div id="71e43f32" class="cell" data-execution_count="1">
+<div id="b8add02b" class="cell" data-execution_count="1">
 <details class="code-fold">
 <summary>Code</summary>
 <div class="sourceCode cell-code" id="cb1"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> numpy <span class="im">as</span> np</span>

reference/node/tabulated-rating-curve.html

@@ -480,7 +480,7 @@ <h2 data-number="1.1" class="anchored" data-anchor-id="static"><span class="head
 <section id="interpolation" class="level3" data-number="1.1.1">
 <h3 data-number="1.1.1" class="anchored" data-anchor-id="interpolation"><span class="header-section-number">1.1.1</span> Interpolation</h3>
 <p>The <span class="math inline">\(Q(h)\)</span> relationship of a tabulated rating curve is defined as a linear interpolation.</p>
-<div id="3330e0b1" class="cell" data-execution_count="1">
+<div id="54414e21" class="cell" data-execution_count="1">
 <details class="code-fold">
 <summary>Code</summary>
 <div class="sourceCode cell-code" id="cb1"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> numpy <span class="im">as</span> np</span>

reference/test-models.html

@@ -343,7 +343,7 @@ <h1 class="title">Test models</h1>
 
 
 <p>Ribasim developers use the following models in their testbench and in order to test new features.</p>
-<div id="e62e1704" class="cell" data-execution_count="1">
+<div id="a5802d03" class="cell" data-execution_count="1">
 <details class="code-fold">
 <summary>Code</summary>
 <div class="sourceCode cell-code" id="cb1"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> ribasim_testmodels</span>

reference/usage.html

@@ -459,8 +459,8 @@ <h1 data-number="1"><span class="header-section-number">1</span> Configuration f
 <span id="cb1-33"><a href="#cb1-33" aria-hidden="true" tabindex="-1"></a><span class="dt">force_dtmin</span> <span class="op">=</span> <span class="cn">false</span> <span class="co"># optional, default false</span></span>
 <span id="cb1-34"><a href="#cb1-34" aria-hidden="true" tabindex="-1"></a><span class="dt">abstol</span> <span class="op">=</span> <span class="fl">1e-6</span>       <span class="co"># optional, default 1e-6</span></span>
 <span id="cb1-35"><a href="#cb1-35" aria-hidden="true" tabindex="-1"></a><span class="dt">reltol</span> <span class="op">=</span> <span class="fl">1e-5</span>       <span class="co"># optional, default 1e-5</span></span>
-<span id="cb1-36"><a href="#cb1-36" aria-hidden="true" tabindex="-1"></a><span class="dt">water_balance_abstol</span> <span class="op">=</span> <span class="fl">1e-6</span> <span class="co"># optional, default 1e-6</span></span>
-<span id="cb1-37"><a href="#cb1-37" aria-hidden="true" tabindex="-1"></a><span class="dt">water_balance_reltol</span> <span class="op">=</span> <span class="fl">1e-6</span> <span class="co"># optional, default 1e-6</span></span>
+<span id="cb1-36"><a href="#cb1-36" aria-hidden="true" tabindex="-1"></a><span class="dt">water_balance_abstol</span> <span class="op">=</span> <span class="fl">1e-3</span> <span class="co"># optional, default 1e-3</span></span>
+<span id="cb1-37"><a href="#cb1-37" aria-hidden="true" tabindex="-1"></a><span class="dt">water_balance_reltol</span> <span class="op">=</span> <span class="fl">1e-2</span> <span class="co"># optional, default 1e-2</span></span>
 <span id="cb1-38"><a href="#cb1-38" aria-hidden="true" tabindex="-1"></a><span class="dt">maxiters</span> <span class="op">=</span> <span class="fl">1e9</span>      <span class="co"># optional, default 1e9</span></span>
 <span id="cb1-39"><a href="#cb1-39" aria-hidden="true" tabindex="-1"></a><span class="dt">sparse</span> <span class="op">=</span> <span class="cn">true</span>       <span class="co"># optional, default true</span></span>
 <span id="cb1-40"><a href="#cb1-40" aria-hidden="true" tabindex="-1"></a><span class="dt">autodiff</span> <span class="op">=</span> <span class="cn">true</span>     <span class="co"># optional, default true</span></span>
@@ -713,6 +713,7 @@ <h2 data-number="5.1" class="anchored" data-anchor-id="basin---basin.arrow"><spa
 <li>The <code>balance_error</code> is the difference between the <code>storage_rate</code> on one side and the <code>inflow_rate</code> and <code>outflow_rate</code> on the other side: <code>storage_rate</code> - (<code>inflow_rate</code> - <code>outflow_rate</code>). It can be used to check if the numerical error when solving the water balance is sufficiently small.</li>
 <li>The <code>relative_error</code> is the fraction of the <code>balance_error</code> over the mean of the <code>total_inflow</code> and <code>total_outflow</code>.</li>
 </ul>
+<p>For a more in-depth explanation of the water balance error see <a href="../concept/equations.html#the-water-balance-error">here</a>.</p>
 <table class="caption-top">
 <thead>
 <tr>