Only gravity driven Outlet flow (#549)

Fixes #548.

Currently the outlet switches on and off in a discrete way due to a
switch in sign of the level difference. This can result in the outlet
rapidly turning on and off:


![flipperende_outlet](https://github.com/Deltares/Ribasim/assets/74617371/d834ffa9-ac48-4fb1-99c5-732369f43c35)

To solve this problem, a reduction factor should be applied to the
outlet flow, based on the level difference. This reduction factor will
have to be taken into account in the implementation of the analytical
Jacobian (with an unsatisfying effort to improvement ratio). This I
think will be a recurring problem until we are able to switch to an AD
Jacobian.

---------

Co-authored-by: Martijn Visser <[email protected]>

build/libribasim/src/libribasim.jl

@@ -27,7 +27,7 @@ This expands to
 ```
 try
     global model
-    isnothing(model) && error("Model not initialized")
+    model === nothing && error("Model not initialized")
     BMI.update(model)
 catch
     Base.invokelatest(Base.display_error, current_exceptions())
@@ -40,7 +40,7 @@ macro try_c(ex)
     return quote
         try
             global model
-            isnothing(model) && error("Model not initialized")
+            model === nothing && error("Model not initialized")
             $(esc(ex))
         catch e
             global last_error_message
@@ -81,7 +81,7 @@ end
 
 Base.@ccallable function finalize()::Cint
     @try_c_uninitialized begin
-        if !isnothing(model)
+        if model !== nothing
             BMI.finalize(model)
         end
         model = nothing

core/src/bmi.jl

@@ -53,7 +53,7 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model
         for id in pid_control.node_id
             id_controlled = only(outneighbors(connectivity.graph_control, id))
             pump_idx = findsorted(pump.node_id, id_controlled)
-            if isnothing(pump_idx)
+            if pump_idx === nothing
                 outlet_idx = findsorted(outlet.node_id, id_controlled)
                 outlet.is_pid_controlled[outlet_idx] = true
             else
@@ -247,7 +247,7 @@ function get_value(
 
         if hasindex_basin
             _, level = get_area_and_level(basin, basin_idx, u[basin_idx])
-        elseif !isnothing(level_boundary_idx)
+        elseif level_boundary_idx !== nothing
             level = level_boundary.level[level_boundary_idx](t + Δt)
         else
             error(
@@ -260,7 +260,7 @@ function get_value(
     elseif variable == "flow_rate"
         flow_boundary_idx = findsorted(flow_boundary.node_id, feature_id)
 
-        if isnothing(flow_boundary_idx)
+        if flow_boundary_idx === nothing
             error("Flow condition node #$feature_id is not a flow boundary.")
         end
 

core/src/config.jl

@@ -84,7 +84,7 @@ end
 
 function Base.convert(::Type{Compression}, str::AbstractString)
     i = findfirst(==(Symbol(str)) ∘ Symbol, instances(Compression))
-    if isnothing(i)
+    if i === nothing
         throw(
             ArgumentError(
                 "Compression algorithm $str not supported, choose one of: $(join(instances(Compression), " ")).",
@@ -149,15 +149,15 @@ function Base.show(io::IO, c::Config)
     println(io, "Ribasim Config")
     for field in fieldnames(typeof(c))
         f = getfield(c, field)
-        isnothing(f) || println(io, "\t$field\t= $f")
+        f === nothing || println(io, "\t$field\t= $f")
     end
 end
 
 function Base.show(io::IO, c::TableOption)
     first = true
     for field in fieldnames(typeof(c))
         f = getfield(c, field)
-        if !isnothing(f)
+        if f !== nothing
             first && (first = false; println(io))
             println(io, "\t\t$field\t= $f")
         end
@@ -180,7 +180,7 @@ const algorithms = Dict{String, Type}(
 "Create an OrdinaryDiffEqAlgorithm from solver config"
 function algorithm(solver::Solver)::OrdinaryDiffEqAlgorithm
     algotype = get(algorithms, solver.algorithm, nothing)
-    if isnothing(algotype)
+    if algotype === nothing
         options = join(keys(algorithms), ", ")
         error("Given solver algorithm $(solver.algorithm) not supported.\n\
             Available options are: ($(options)).")

core/src/create.jl

@@ -74,14 +74,14 @@ function parse_static_and_time(
     end
 
     # Get node IDs of static nodes if the static table exists
-    static_node_ids = if isnothing(static)
+    static_node_ids = if static === nothing
         Set{Int}()
     else
         Set(static.node_id)
     end
 
     # Get node IDs of transient nodes if the time table exists
-    time_node_ids = if isnothing(time)
+    time_node_ids = if time === nothing
         Set{Int}()
     else
         Set(time.node_id)
@@ -433,7 +433,8 @@ end
 
 function Outlet(db::DB, config::Config, chunk_size::Int)::Outlet
     static = load_structvector(db, config, OutletStaticV1)
-    defaults = (; min_flow_rate = 0.0, max_flow_rate = NaN, active = true)
+    defaults =
+        (; min_flow_rate = 0.0, max_flow_rate = NaN, min_crest_level = NaN, active = true)
     parsed_parameters, valid = parse_static_and_time(db, config, "Outlet"; static, defaults)
     is_pid_controlled = falses(length(parsed_parameters.node_id))
 
@@ -454,6 +455,7 @@ function Outlet(db::DB, config::Config, chunk_size::Int)::Outlet
         flow_rate,
         parsed_parameters.min_flow_rate,
         parsed_parameters.max_flow_rate,
+        parsed_parameters.min_crest_level,
         parsed_parameters.control_mapping,
         is_pid_controlled,
     )

core/src/io.jl

@@ -50,7 +50,7 @@ function load_data(
     path = getfield(getfield(config, snake_case(node)), kind)
     sqltable = tablename(schema)
 
-    table = if !isnothing(path)
+    table = if path !== nothing
         table_path = input_path(config, path)
         Table(read(table_path))
     elseif exists(db, sqltable)
@@ -76,7 +76,7 @@ function load_structvector(
 )::StructVector{T} where {T <: AbstractRow}
     table = load_data(db, config, T)
 
-    if isnothing(table)
+    if table === nothing
         return StructVector{T}(undef, 0)
     end
 
@@ -98,7 +98,7 @@ function load_structvector(
     table = StructVector{T}(nt)
     sv = Legolas._schema_version_from_record_type(T)
     tableschema = Tables.schema(table)
-    if declared(sv) && !isnothing(tableschema)
+    if declared(sv) && tableschema !== nothing
         validate(tableschema, sv)
         # R = Legolas.record_type(sv)
         # foreach(R, Tables.rows(table))  # construct each row