From 03c57b0767a7d6cba9580c1dd8b6b714d7ffc9d9 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Wed, 4 Dec 2024 13:30:54 +0800 Subject: [PATCH] move geometry to parameters --- src/erosion.jl | 4 +- src/erosion/overland_flow_erosion.jl | 2 +- src/erosion/rainfall_erosion.jl | 4 +- src/erosion/river_erosion.jl | 2 +- src/geometry.jl | 79 ------------------- src/parameters.jl | 82 ++++++++++++++++++++ src/sediment_flux.jl | 8 +- src/sediment_transport/river_transport.jl | 4 +- src/sediment_transport/transport_capacity.jl | 12 +-- test/sediment_config.toml | 4 +- 10 files changed, 102 insertions(+), 99 deletions(-) delete mode 100644 src/geometry.jl diff --git a/src/erosion.jl b/src/erosion.jl index c6951587c..c3721132c 100644 --- a/src/erosion.jl +++ b/src/erosion.jl @@ -1,6 +1,6 @@ @with_kw struct SoilLoss{RE, OFE, SE, T} hydrometeo_forcing::HydrometeoForcing - geometry::LandGeometry + geometry::LandParameters rainfall_erosion::RE overland_flow_erosion::OFE soil_erosion::SE @@ -10,7 +10,7 @@ function SoilLoss(nc, config, inds) n = length(inds) hydrometeo_forcing = HydrometeoForcing(n) - geometry = LandGeometry(nc, config, inds) + geometry = LandParameters(nc, config, inds) # Rainfall erosion rainfallerosionmodel = get(config.model, "rainfall_erosion", "answers")::String diff --git a/src/erosion/overland_flow_erosion.jl b/src/erosion/overland_flow_erosion.jl index 494dc2ff2..1b82968d5 100644 --- a/src/erosion/overland_flow_erosion.jl +++ b/src/erosion/overland_flow_erosion.jl @@ -99,7 +99,7 @@ function update_boundary_conditions!( return nothing end -function update!(model::OverlandFlowErosionAnswersModel, geometry::LandGeometry, ts) +function update!(model::OverlandFlowErosionAnswersModel, geometry::LandParameters, ts) (; q) = model.boundary_conditions (; usle_k, usle_c, answers_k) = model.parameters (; amount) = model.variables diff --git a/src/erosion/rainfall_erosion.jl b/src/erosion/rainfall_erosion.jl index 5e70a5aaa..66ec3479c 100644 --- a/src/erosion/rainfall_erosion.jl +++ b/src/erosion/rainfall_erosion.jl @@ -127,7 +127,7 @@ function update_boundary_conditions!( @. waterlevel = model.boundary_conditions.waterlevel_land end -function update!(model::RainfallErosionEurosemModel, geometry::LandGeometry, ts) +function update!(model::RainfallErosionEurosemModel, geometry::LandParameters, ts) (; precipitation, interception, waterlevel) = model.boundary_conditions (; soil_detachability, @@ -221,7 +221,7 @@ function update_boundary_conditions!( @. precipitation = hydrometeo_forcing.precipitation end -function update!(model::RainfallErosionAnswersModel, geometry::LandGeometry, ts) +function update!(model::RainfallErosionAnswersModel, geometry::LandParameters, ts) (; precipitation) = model.boundary_conditions (; usle_k, usle_c) = model.parameters (; amount) = model.variables diff --git a/src/erosion/river_erosion.jl b/src/erosion/river_erosion.jl index b6f99ffcf..655bed9b2 100644 --- a/src/erosion/river_erosion.jl +++ b/src/erosion/river_erosion.jl @@ -73,7 +73,7 @@ function update_boundary_conditions!( @. waterlevel = waterlevel_river end -function update!(model::RiverErosionJulianTorresModel, geometry::RiverGeometry, ts) +function update!(model::RiverErosionJulianTorresModel, geometry::RiverParameters, ts) (; waterlevel) = model.boundary_conditions (; d50) = model.parameters (; bed, bank) = model.variables diff --git a/src/geometry.jl b/src/geometry.jl deleted file mode 100644 index 476731a7e..000000000 --- a/src/geometry.jl +++ /dev/null @@ -1,79 +0,0 @@ - -@get_units @grid_loc @with_kw struct LandGeometry{T} - # cell area [m^2] - area::Vector{T} | "m^2" - # drain width [m] - width::Vector{T} | "m" - # drain slope - slope::Vector{T} | "-" -end - -function LandGeometry(nc, config, inds) - # read x, y coordinates and calculate cell length [m] - y_nc = read_y_axis(nc) - x_nc = read_x_axis(nc) - y = permutedims(repeat(y_nc; outer = (1, length(x_nc))))[inds] - cellength = abs(mean(diff(x_nc))) - - sizeinmetres = get(config.model, "sizeinmetres", false)::Bool - xl, yl = cell_lengths(y, cellength, sizeinmetres) - area = xl .* yl - ldd = ncread(nc, config, "ldd"; optional = false, sel = inds, allow_missing = true) - drain_width = map(detdrainwidth, ldd, xl, yl) - landslope = ncread( - nc, - config, - "vertical.geometry.slope"; - optional = false, - sel = inds, - type = Float, - ) - clamp!(landslope, 0.00001, Inf) - - land_geometry = - LandGeometry{Float}(; area = area, width = drain_width, slope = landslope) - return land_geometry -end - -@get_units @grid_loc @with_kw struct RiverGeometry{T} - # drain width [m] - width::Vector{T} | "m" - # drain length - length::Vector{T} | "m" - # slope - slope::Vector{T} | "-" -end - -function RiverGeometry(nc, config, inds) - riverwidth = ncread( - nc, - config, - "lateral.river.geometry.width"; - optional = false, - sel = inds, - type = Float, - ) - riverlength = ncread( - nc, - config, - "lateral.river.geometry.length"; - optional = false, - sel = inds, - type = Float, - ) - riverslope = ncread( - nc, - config, - "lateral.river.geometry.slope"; - optional = false, - sel = inds, - type = Float, - ) - minimum(riverlength) > 0 || error("river length must be positive on river cells") - minimum(riverwidth) > 0 || error("river width must be positive on river cells") - clamp!(riverslope, 0.00001, Inf) - - river_geometry = - RiverGeometry{Float}(; width = riverwidth, length = riverlength, slope = riverslope) - return river_geometry -end \ No newline at end of file diff --git a/src/parameters.jl b/src/parameters.jl index f2d2d582a..b6f1d085e 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -101,4 +101,86 @@ function VegetationParameters(nc, config, inds) ) end return vegetation_parameter_set +end + +@get_units @grid_loc @with_kw struct LandParameters{T} + # cell area [m^2] + area::Vector{T} | "m^2" + # drain width [m] + width::Vector{T} | "m" + # drain slope + slope::Vector{T} | "-" +end + +function LandParameters(nc, config, inds) + # read x, y coordinates and calculate cell length [m] + y_nc = read_y_axis(nc) + x_nc = read_x_axis(nc) + y = permutedims(repeat(y_nc; outer = (1, length(x_nc))))[inds] + cellength = abs(mean(diff(x_nc))) + + sizeinmetres = get(config.model, "sizeinmetres", false)::Bool + xl, yl = cell_lengths(y, cellength, sizeinmetres) + area = xl .* yl + ldd = ncread(nc, config, "ldd"; optional = false, sel = inds, allow_missing = true) + drain_width = map(detdrainwidth, ldd, xl, yl) + landslope = ncread( + nc, + config, + "vertical.land_parameter_set.slope"; + optional = false, + sel = inds, + type = Float, + ) + clamp!(landslope, 0.00001, Inf) + + land_parameter_set = + LandParameters{Float}(; area = area, width = drain_width, slope = landslope) + return land_parameter_set +end + +@get_units @grid_loc @with_kw struct RiverParameters{T} + # river width [m] + width::Vector{T} | "m" + # river length + length::Vector{T} | "m" + # slope + slope::Vector{T} | "-" +end + +function RiverParameters(nc, config, inds) + riverwidth = ncread( + nc, + config, + "lateral.river_parameter_set.width"; + optional = false, + sel = inds, + type = Float, + ) + riverlength = ncread( + nc, + config, + "lateral.river_parameter_set.length"; + optional = false, + sel = inds, + type = Float, + ) + riverslope = ncread( + nc, + config, + "lateral.river_parameter_set.slope"; + optional = false, + sel = inds, + type = Float, + ) + minimum(riverlength) > 0 || error("river length must be positive on river cells") + minimum(riverwidth) > 0 || error("river width must be positive on river cells") + clamp!(riverslope, 0.00001, Inf) + + river_parameter_set = RiverParameters{Float}(; + width = riverwidth, + length = riverlength, + slope = riverslope, + ) + return river_parameter_set end \ No newline at end of file diff --git a/src/sediment_flux.jl b/src/sediment_flux.jl index bcb710170..2560ebd23 100644 --- a/src/sediment_flux.jl +++ b/src/sediment_flux.jl @@ -1,7 +1,7 @@ ### Overland flow ### @get_units @grid_loc @with_kw struct OverlandFlowSediment{TT, SF, TR, T} hydrometeo_forcing::HydrometeoForcing - geometry::LandGeometry + geometry::LandParameters transport_capacity::TT sediment_flux::SF to_river::TR @@ -12,7 +12,7 @@ end function OverlandFlowSediment(nc, config, inds, waterbodies, rivers) n = length(inds) hydrometeo_forcing = HydrometeoForcing(n) - geometry = LandGeometry(nc, config, inds) + geometry = LandParameters(nc, config, inds) # Check what transport capacity equation will be used do_river = get(config.model, "runrivermodel", false)::Bool # Overland flow transport capacity method: ["yalinpart", "govers", "yalin"] @@ -80,7 +80,7 @@ end ### River ### @get_units @grid_loc @with_kw struct RiverSediment{TTR, ER, SFR, CR, T} hydrometeo_forcing::HydrometeoForcing - geometry::RiverGeometry + geometry::RiverParameters transport_capacity::TTR potential_erosion::ER sediment_flux::SFR @@ -91,7 +91,7 @@ end function RiverSediment(nc, config, inds, waterbodies) n = length(inds) hydrometeo_forcing = HydrometeoForcing(n) - geometry = RiverGeometry(nc, config, inds) + geometry = RiverParameters(nc, config, inds) # Check what transport capacity equation will be used # River flow transport capacity method: ["bagnold", "engelund", "yang", "kodatie", "molinas"] diff --git a/src/sediment_transport/river_transport.jl b/src/sediment_transport/river_transport.jl index 867ebbc38..73abfc9a1 100644 --- a/src/sediment_transport/river_transport.jl +++ b/src/sediment_transport/river_transport.jl @@ -385,7 +385,7 @@ end function update!( model::SedimentRiverTransportModel, network, - geometry::RiverGeometry, + geometry::RiverParameters, waterbodies, ts, ) @@ -952,7 +952,7 @@ function update_boundary_conditions!( @. gravel = sediment_flux_model.variables.gravel end -function update!(model::SedimentConcentrationsRiverModel, geometry::RiverGeometry, ts) +function update!(model::SedimentConcentrationsRiverModel, geometry::RiverParameters, ts) (; q, waterlevel, clay, silt, sand, sagg, lagg, gravel) = model.boundary_conditions (; dm_clay, dm_silt, dm_sand, dm_sagg, dm_lagg, dm_gravel) = model.parameters (; total, suspended, bed) = model.variables diff --git a/src/sediment_transport/transport_capacity.jl b/src/sediment_transport/transport_capacity.jl index e1f625dfa..e5fdfb228 100644 --- a/src/sediment_transport/transport_capacity.jl +++ b/src/sediment_transport/transport_capacity.jl @@ -355,7 +355,7 @@ end function update!( model::TransportCapacityYalinDifferentiationModel, - geometry::LandGeometry, + geometry::LandParameters, waterbodies, rivers, ts, @@ -520,7 +520,7 @@ function TransportCapacityBagnoldModel(nc, config, inds) return model end -function update!(model::TransportCapacityBagnoldModel, geometry::RiverGeometry, ts) +function update!(model::TransportCapacityBagnoldModel, geometry::RiverParameters, ts) (; q, waterlevel) = model.boundary_conditions (; c_bagnold, e_bagnold) = model.parameters (; amount) = model.variables @@ -561,7 +561,7 @@ function TransportCapacityEngelundModel(nc, config, inds) return model end -function update!(model::TransportCapacityEngelundModel, geometry::RiverGeometry, ts) +function update!(model::TransportCapacityEngelundModel, geometry::RiverParameters, ts) (; q, waterlevel) = model.boundary_conditions (; density, d50) = model.parameters (; amount) = model.variables @@ -655,7 +655,7 @@ function TransportCapacityKodatieModel(nc, config, inds) return model end -function update!(model::TransportCapacityKodatieModel, geometry::RiverGeometry, ts) +function update!(model::TransportCapacityKodatieModel, geometry::RiverParameters, ts) (; q, waterlevel) = model.boundary_conditions (; a_kodatie, b_kodatie, c_kodatie, d_kodatie) = model.parameters (; amount) = model.variables @@ -697,7 +697,7 @@ function TransportCapacityYangModel(nc, config, inds) return model end -function update!(model::TransportCapacityYangModel, geometry::RiverGeometry, ts) +function update!(model::TransportCapacityYangModel, geometry::RiverParameters, ts) (; q, waterlevel) = model.boundary_conditions (; density, d50) = model.parameters (; amount) = model.variables @@ -737,7 +737,7 @@ function TransportCapacityMolinasModel(nc, config, inds) return model end -function update!(model::TransportCapacityMolinasModel, geometry::RiverGeometry, ts) +function update!(model::TransportCapacityMolinasModel, geometry::RiverParameters, ts) (; q, waterlevel) = model.boundary_conditions (; density, d50) = model.parameters (; amount) = model.variables diff --git a/test/sediment_config.toml b/test/sediment_config.toml index de64a3188..0efaf68ab 100644 --- a/test/sediment_config.toml +++ b/test/sediment_config.toml @@ -79,7 +79,7 @@ waterlevel_land = "levKinL" #interception = "int" q_land = "runL" -[input.vertical.geometry] +[input.vertical.land_parameter_set] slope = "Slope" [input.vertical.rainfall_erosion.parameters] @@ -122,7 +122,7 @@ dm_lagg = "dm_lagg" waterlevel_river = "h" q_river = "q" -[input.lateral.river.geometry] +[input.lateral.river_parameter_set] length = "wflow_riverlength" slope = "RiverSlope" width = "wflow_riverwidth"