Repository: Deltares/Wflow.jl (julia)

Issue:
BMI and request of scalar variables
### Wflow version checks

- [X] I have checked that this issue has not already been reported.

- [X] I have checked that this bug exists on the latest version of Wflow.


### Reproducible Example

`BMI.get_var_itemsize(model, "lateral.land.alpha_pow")`

### Current behaviour

```
{"status":"ERROR","error":"Wflow BMI get_var_itemsize failed
exception =
MethodError: no method matching view(::Float64, ::UnitRange{Int64}
```

### Desired behaviour

Should give the item size (8) back. The error is thrown because a `view` is applied on a Float64 in `BMI.get_value_ptr(model, name)`. For now, decided to not exchange scalar variables through BMI (there are only a few model parameters defined on a scalar grid).

### Additional Context

_No response_

Pull request description:
## Issue addressed
Fixes #331 

## Explanation
See issue and updated changelog.

## Checklist
- [x] Updated tests or added new tests
- [x] Branch is up to date with `master`
- [x] Tests & pre-commit hooks pass
- [x] Updated documentation if needed
- [x] Updated changelog.md if needed

Patch:
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index c14da9a0c..93abfcc4d 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -18,6 +18,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   per grid cell in millimeters. Excluded from this variable are the floodplain, lakes and
   reservoirs.
 
+### Changed
+- Stop exposing scalar variables through BMI. The `BMI.get_value_ptr` function was
+  not working correctly for scalar model variables (a `view` was applied). Only a few scalar
+  model parameters are defined, and it is not expected that exposing these variables is
+  required (e.g. for model coupling) while code changes for these variables (including
+  struct fields) are required.
+
 ## v0.7.3 - 2024-01-12
 
 ### Fixed
diff --git a/src/bmi.jl b/src/bmi.jl
index 0ef6e574a..1c7d8f476 100644
--- a/src/bmi.jl
+++ b/src/bmi.jl
@@ -4,13 +4,13 @@
 # Mapping of grid identifier to a key, to get the active indices of the model domain.
 # See also function active_indices(network, key::Tuple).
 const grids = Dict{Int,Tuple{Symbol}}(
-    1 => (:reservoir,),
-    2 => (:lake,),
-    3 => (:drain,),
-    4 => (:river,),
+    0 => (:reservoir,),
+    1 => (:lake,),
+    2 => (:drain,),
+    3 => (:river,),
+    4 => (:land,),
     5 => (:land,),
     6 => (:land,),
-    7 => (:land,),
 )
 
 """
@@ -145,22 +145,20 @@ function BMI.get_var_grid(model::Model, name::String)
     if exchange(param(model, key[1:end-1]), key[end]) == 1
         gridtype = grid_type(param(model, key))
         type = typeof(param(model, key[1:end-1]))
-        return if gridtype == "scalar"
+        return if :reservoir in key
             0
-        elseif :reservoir in key
-            1
         elseif :lake in key
-            2
+            1
         elseif :drain in key
-            3
+            2
         elseif :river in key
-            4
+            3
         elseif type <: ShallowWaterLand && occursin("x", s[end])
-            5
+            4
         elseif type <: ShallowWaterLand && occursin("y", s[end])
-            6
+            5
         else
-            7
+            6
         end
     else
         error("$name not listed as variable for BMI exchange")
@@ -293,11 +291,9 @@ function BMI.set_value_at_indices(
 end
 
 function BMI.get_grid_type(model::Model, grid::Int)
-    if grid == 0
-        "scalar"
-    elseif grid in 1:3
+    if grid in 0:2
         "points"
-    elseif grid in 4:7
+    elseif grid in 3:6
         "unstructured"
     else
         error("unknown grid type $grid")
@@ -305,9 +301,7 @@ function BMI.get_grid_type(model::Model, grid::Int)
 end
 
 function BMI.get_grid_rank(model::Model, grid::Int)
-    if grid == 0
-        0
-    elseif grid in 1:7
+    if grid in 0:6
         2
     else
         error("unknown grid type $grid")
@@ -344,13 +338,13 @@ end
 
 function BMI.get_grid_edge_count(model::Model, grid::Int)
     @unpack network = model
-    if grid == 4
+    if grid == 3
         return ne(network.river.graph)
-    elseif grid == 5
+    elseif grid == 4
         return length(network.land.staggered_indices.xu)
-    elseif grid == 6
+    elseif grid == 5
         return length(network.land.staggered_indices.yu)
-    elseif grid in 0:3 || grid == 7
+    elseif grid in 0:2 || grid == 6
         warn("edges are not provided for grid type $grid (variables are located at nodes)")
     else
         error("unknown grid type $grid")
@@ -362,25 +356,25 @@ function BMI.get_grid_edge_nodes(model::Model, grid::Int, edge_nodes::Vector{Int
     n = length(edge_nodes)
     m = div(n, 2)
     # inactive nodes (boundary/ghost points) are set at -999
-    if grid == 4
+    if grid == 3
         nodes_at_edge = adjacent_nodes_at_link(network.river.graph)
         nodes_at_edge.dst[nodes_at_edge.dst.==m+1] .= -999
         edge_nodes[range(1, n, step = 2)] = nodes_at_edge.src
         edge_nodes[range(2, n, step = 2)] = nodes_at_edge.dst
         return edge_nodes
-    elseif grid == 5
+    elseif grid == 4
         xu = network.land.staggered_indices.xu
         edge_nodes[range(1, n, step = 2)] = 1:m
         xu[xu.==m+1] .= -999
         edge_nodes[range(2, n, step = 2)] = xu
         return edge_nodes
-    elseif grid == 6
+    elseif grid == 5
         yu = network.land.staggered_indices.yu
         edge_nodes[range(1, n, step = 2)] = 1:m
         yu[yu.==m+1] .= -999
         edge_nodes[range(2, n, step = 2)] = yu
         return edge_nodes
-    elseif grid in 0:3 || grid == 7
+    elseif grid in 0:2 || grid == 6
         @warn("edges are not provided for grid type $grid (variables are located at nodes)")
     else
         error("unknown grid type $grid")
diff --git a/src/flow.jl b/src/flow.jl
index af77daa89..c5e81ad5d 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -2,7 +2,7 @@
 abstract type SurfaceFlow end
 
 @get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowRiver{T,R,L} <: SurfaceFlow
-    β::T | "-" | _ | "scalar"                    # constant in Manning's equation
+    β::T | "-" | 0 | "scalar"                    # constant in Manning's equation
     sl::Vector{T} | "m m-1"                      # Slope [m m⁻¹]
     n::Vector{T} | "s m-1/3"                     # Manning's roughness [s m⁻⅓]
     dl::Vector{T} | "m"                          # Drain length [m]
@@ -19,7 +19,7 @@ abstract type SurfaceFlow end
     Δt::T | "s" | 0 | "none" | "none"            # Model time step [s]
     its::Int | "-" | 0 | "none" | "none"         # Number of fixed iterations
     width::Vector{T} | "m"                       # Flow width [m]
-    alpha_pow::T | "-" | _ | "scalar"            # Used in the power part of α
+    alpha_pow::T | "-" | 0 | "scalar"            # Used in the power part of α
     alpha_term::Vector{T} | "-"                  # Term used in computation of α
     α::Vector{T} | "s3/5 m1/5"                   # Constant in momentum equation A = αQᵝ, based on Manning's equation
     cel::Vector{T} | "m s-1"                     # Celerity of the kinematic wave
@@ -37,7 +37,7 @@ abstract type SurfaceFlow end
 end
 
 @get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowLand{T} <: SurfaceFlow
-    β::T | "-"  | _ | "scalar"                      # constant in Manning's equation
+    β::T | "-"  | 0 | "scalar"                      # constant in Manning's equation
     sl::Vector{T} | "m m-1"                         # Slope [m m⁻¹]
     n::Vector{T} | "s m-1/3"                        # Manning's roughness [s m⁻⅓]
     dl::Vector{T} | "m"                             # Drain length [m]
@@ -52,7 +52,7 @@ end
     Δt::T | "s" | 0 | "none" | "none"               # Model time step [s]
     its::Int | "-" | 0 | "none" | "none"            # Number of fixed iterations
     width::Vector{T} | "m"                          # Flow width [m]
-    alpha_pow::T | "-" | _ | "scalar"               # Used in the power part of α
+    alpha_pow::T | "-" | 0 | "scalar"               # Used in the power part of α
     alpha_term::Vector{T} | "-"                     # Term used in computation of α
     α::Vector{T} | "s3/5 m1/5"                      # Constant in momentum equation A = αQᵝ, based on Manning's equation
     cel::Vector{T} | "m s-1"                        # Celerity of the kinematic wave
@@ -468,9 +468,9 @@ end
     ne::Int | "-" | 0 | "none" | "none"                     # number of edges/links
     active_n::Vector{Int} | "-"                             # active nodes
     active_e::Vector{Int} | "-" | _ | "edge"                # active edges/links
-    g::T | "m s-2" | _ | "scalar"                           # acceleration due to gravity
-    α::T | "-" | _ | "scalar"                               # stability coefficient (Bates et al., 2010)
-    h_thresh::T | "m" | _ | "scalar"                        # depth threshold for calculating flow
+    g::T | "m s-2" | 0 | "scalar"                           # acceleration due to gravity
+    α::T | "-" | 0 | "scalar"                               # stability coefficient (Bates et al., 2010)
+    h_thresh::T | "m" | 0 | "scalar"                        # depth threshold for calculating flow
     Δt::T | "s" | 0 | "none" | "none"                       # model time step [s]
     q::Vector{T} | "m3 s-1" | _ | "edge"                    # river discharge (subgrid channel)
     q0::Vector{T} | "m3 s-1"| _ | "edge"                    # river discharge (subgrid channel) at previous time step
@@ -951,11 +951,11 @@ const dirs = (:yd, :xd, :xu, :yu)
     yl::Vector{T} | "m"                                     # cell length y direction
     xwidth::Vector{T} | "m" | _ | "edge"                    # effective flow width x direction (floodplain)
     ywidth::Vector{T} | "m" | _ | "edge"                    # effective flow width y direction (floodplain)
-    g::T | "m2 s-1" | _ | "scalar"                          # acceleration due to gravity
-    θ::T | "-" | _ | "scalar"                               # weighting factor (de Almeida et al., 2012)
-    α::T | "-" | _ | "scalar"                               # stability coefficient (de Almeida et al., 2012)
-    h_thresh::T | "m" | _ | "scalar"                        # depth threshold for calculating flow
-    Δt::T | "s" | _ | "scalar"                              # model time step [s]
+    g::T | "m2 s-1" | 0 | "scalar"                          # acceleration due to gravity
+    θ::T | "-" | 0 | "scalar"                               # weighting factor (de Almeida et al., 2012)
+    α::T | "-" | 0 | "scalar"                               # stability coefficient (de Almeida et al., 2012)
+    h_thresh::T | "m" | 0 | "scalar"                        # depth threshold for calculating flow
+    Δt::T | "s" | 0 | "none" | "none"                       # model time step [s]
     qy0::Vector{T} | "m3 s-1" | _ | "edge"                  # flow in y direction at previous time step
     qx0::Vector{T} | "m3 s-1" | _ | "edge"                  # flow in x direction at previous time step
     qx::Vector{T} | "m3 s-1" | _ | "edge"                   # flow in x direction


Test patch:
diff --git a/test/bmi.jl b/test/bmi.jl
index 36154e1b6..13713cd1b 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -20,20 +20,20 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 
         @testset "model information functions" begin
             @test BMI.get_component_name(model) == "sbm"
-            @test BMI.get_input_item_count(model) == 184
-            @test BMI.get_output_item_count(model) == 184
+            @test BMI.get_input_item_count(model) == 180
+            @test BMI.get_output_item_count(model) == 180
             @test BMI.get_input_var_names(model)[[1, 5, 151, 175]] == [
                 "vertical.nlayers",
                 "vertical.θᵣ",
-                "lateral.river.sl",
-                "lateral.river.reservoir.targetminfrac",
+                "lateral.river.q",
+                "lateral.river.reservoir.outflow",
             ]
         end
 
         @testset "variable information functions" begin
-            @test BMI.get_var_grid(model, "vertical.θₛ") == 7
-            @test BMI.get_var_grid(model, "lateral.river.h") == 4
-            @test BMI.get_var_grid(model, "lateral.river.reservoir.inflow") == 1
+            @test BMI.get_var_grid(model, "vertical.θₛ") == 6
+            @test BMI.get_var_grid(model, "lateral.river.h") == 3
+            @test BMI.get_var_grid(model, "lateral.river.reservoir.inflow") == 0
             @test_throws ErrorException BMI.get_var_grid(model, "lateral.river.lake.volume")
             @test BMI.get_var_type(model, "lateral.river.reservoir.inflow") == "$Float"
             @test BMI.get_var_units(model, "vertical.θₛ") == "-"
@@ -41,6 +41,9 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             @test BMI.get_var_nbytes(model, "lateral.river.q") ==
                   length(model.lateral.river.q) * sizeof(Float)
             @test BMI.get_var_location(model, "lateral.river.q") == "node"
+            @test_throws ErrorException(
+                "lateral.land.alpha_pow not listed as variable for BMI exchange",
+            ) BMI.get_var_itemsize(model, "lateral.land.alpha_pow")
         end
 
         model = BMI.update(model)
@@ -79,32 +82,32 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
         end
 
         @testset "model grid functions" begin
-            @test BMI.get_grid_type(model, 0) == "scalar"
+            @test BMI.get_grid_type(model, 0) == "points"
             @test BMI.get_grid_type(model, 2) == "points"
-            @test BMI.get_grid_type(model, 7) == "unstructured"
+            @test BMI.get_grid_type(model, 6) == "unstructured"
             @test_throws ErrorException BMI.get_grid_rank(model, 8)
-            @test BMI.get_grid_rank(model, 0) == 0
-            @test BMI.get_grid_rank(model, 7) == 2
+            @test BMI.get_grid_rank(model, 0) == 2
+            @test BMI.get_grid_rank(model, 6) == 2
             @test_throws ErrorException BMI.get_grid_rank(model, 8)
-            @test BMI.get_grid_node_count(model, 1) == 2
-            @test BMI.get_grid_node_count(model, 4) == 5809
+            @test BMI.get_grid_node_count(model, 0) == 2
+            @test BMI.get_grid_node_count(model, 3) == 5809
+            @test BMI.get_grid_node_count(model, 4) == 50063
             @test BMI.get_grid_node_count(model, 5) == 50063
-            @test BMI.get_grid_node_count(model, 6) == 50063
-            @test BMI.get_grid_size(model, 1) == 2
-            @test BMI.get_grid_size(model, 4) == 5809
+            @test BMI.get_grid_size(model, 0) == 2
+            @test BMI.get_grid_size(model, 3) == 5809
+            @test BMI.get_grid_size(model, 4) == 50063
             @test BMI.get_grid_size(model, 5) == 50063
-            @test BMI.get_grid_size(model, 6) == 50063
             @test minimum(BMI.get_grid_x(model, 5, zeros(Float, 50063))) ≈
                   5.426666666666667f0
             @test maximum(BMI.get_grid_x(model, 5, zeros(Float, 50063))) ≈
                   7.843333333333344f0
-            @test BMI.get_grid_x(model, 1, zeros(Float, 2)) ≈
+            @test BMI.get_grid_x(model, 0, zeros(Float, 2)) ≈
                   [5.760000000000002f0, 5.918333333333336f0]
-            @test BMI.get_grid_y(model, 1, zeros(Float, 2)) ≈
+            @test BMI.get_grid_y(model, 0, zeros(Float, 2)) ≈
                   [48.92583333333333f0, 49.909166666666664f0]
-            @test BMI.get_grid_node_count(model, 1) == 2
-            @test BMI.get_grid_edge_count(model, 4) == 5808
-            @test BMI.get_grid_edge_nodes(model, 4, fill(0, 2 * 5808))[1:6] ==
+            @test BMI.get_grid_node_count(model, 0) == 2
+            @test BMI.get_grid_edge_count(model, 3) == 5808
+            @test BMI.get_grid_edge_nodes(model, 3, fill(0, 2 * 5808))[1:6] ==
                   [1, 5, 2, 1, 3, 2]
         end
 
@@ -120,17 +123,17 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
     @testset "BMI grid edges" begin
         tomlpath = joinpath(@__DIR__, "sbm_swf_config.toml")
         model = BMI.initialize(Wflow.Model, tomlpath)
-        @test BMI.get_var_grid(model, "lateral.land.qx") == 5
-        @test BMI.get_var_grid(model, "lateral.land.qy") == 6
+        @test BMI.get_var_grid(model, "lateral.land.qx") == 4
+        @test BMI.get_var_grid(model, "lateral.land.qy") == 5
+        @test BMI.get_grid_edge_count(model, 4) == 50063
         @test BMI.get_grid_edge_count(model, 5) == 50063
-        @test BMI.get_grid_edge_count(model, 6) == 50063
-        @test BMI.get_grid_edge_nodes(model, 5, fill(0, 2 * 50063))[1:4] == [1, -999, 2, 3]
-        @test BMI.get_grid_edge_nodes(model, 6, fill(0, 2 * 50063))[1:4] == [1, 4, 2, 10]
+        @test BMI.get_grid_edge_nodes(model, 4, fill(0, 2 * 50063))[1:4] == [1, -999, 2, 3]
+        @test BMI.get_grid_edge_nodes(model, 5, fill(0, 2 * 50063))[1:4] == [1, 4, 2, 10]
         @test_logs (
             :warn,
-            "edges are not provided for grid type 3 (variables are located at nodes)",
-        ) BMI.get_grid_edge_nodes(model, 3, fill(0, 2 * 50063))
-        @test_throws ErrorException BMI.get_grid_edge_nodes(model, 8, fill(0, 2 * 50063))
+            "edges are not provided for grid type 2 (variables are located at nodes)",
+        ) BMI.get_grid_edge_nodes(model, 2, fill(0, 2 * 50063))
+        @test_throws ErrorException BMI.get_grid_edge_nodes(model, 7, fill(0, 2 * 50063))
         BMI.finalize(model)
     end