|
| 1 | +export update_params!, optim_params |
| 2 | + |
| 3 | +## |
| 4 | +@with_kw mutable struct SoilColumn{FT<:AbstractFloat,N, |
| 5 | + H<:HydraulicProfile{FT,N},T<:ThermalProfile{FT,N}} |
| 6 | + hydraulic::H # 水力剖面 |
| 7 | + thermal::T # 热力剖面 |
| 8 | +end |
| 9 | + |
| 10 | +# 外构造器：N 作为类型参数显式传入，避免 runtime N → 类型不稳定 |
| 11 | +function SoilColumn{FT,N}( |
| 12 | + hydraulic=HydraulicProfile{FT,N}(), |
| 13 | + thermal=ThermalProfile{FT,N}()) where {FT<:AbstractFloat,N} |
| 14 | + |
| 15 | + SoilColumn{FT,N,typeof(hydraulic),typeof(thermal)}(hydraulic, thermal) |
| 16 | +end |
| 17 | + |
| 18 | +# effective_ksat(ps::SoilColumn, i::Int, z_cm::Real) = kv_at_depth(ps.hydraulic.kv, i, z_cm) |
| 19 | + |
| 20 | +# Sync Ksat from kv profile into the SoA hydraulic profile (writes to profile.Ksat[i]). |
| 21 | +function _sync_ksat!(kv::Union{AbstractKv,AbstractKvLayers}, |
| 22 | + profile::AbstractRetentionLayers{FT,N}, |
| 23 | + dz_cm::AbstractVector) where {FT,N} |
| 24 | + z_cm = FT(0) |
| 25 | + @inbounds for i in 1:N |
| 26 | + dz_i = isempty(dz_cm) ? FT(0) : FT(dz_cm[i]) |
| 27 | + profile.Ksat[i] = kv_layer_ksat(kv, i, z_cm, z_cm + dz_i) |
| 28 | + z_cm += dz_i |
| 29 | + end |
| 30 | +end |
| 31 | + |
| 32 | + |
| 33 | +# mod: :hydraulic | :thermal | [:hydraulic, :thermal] | :all |
| 34 | +# list_fix: field names excluded from calibration, matched via params.name |
| 35 | +# list_sameLayer: field names shared across all hydraulic layers (e.g. Campbell ψ_sat); |
| 36 | +# returns one deduped row per name; update_params! broadcasts it back. |
| 37 | +function optim_params(ps::SoilColumn, mod::Union{Symbol,AbstractVector{Symbol}}; |
| 38 | + inds=nothing, |
| 39 | + list_sameLayer::Vector{Symbol}=Symbol[], |
| 40 | + list_fix::Vector{Symbol}=Symbol[]) |
| 41 | + |
| 42 | + params = parameters(ps) # base struct-traversal in Params.jl |
| 43 | + n = nrow(params) |
| 44 | + mask = trues(n) |
| 45 | + |
| 46 | + # module filter |
| 47 | + if mod !== :all |
| 48 | + mods = mod isa Symbol ? (mod,) : Tuple(mod) |
| 49 | + for (i, p) in enumerate(params.path) |
| 50 | + mask[i] && p[1] ∉ mods && (mask[i] = false) |
| 51 | + end |
| 52 | + end |
| 53 | + |
| 54 | + # layer index filter |
| 55 | + if !isnothing(inds) |
| 56 | + inds_set = Set(inds) |
| 57 | + for (i, p) in enumerate(params.path) |
| 58 | + mask[i] && isa(p[end], Integer) && p[end] ∉ inds_set && (mask[i] = false) |
| 59 | + end |
| 60 | + end |
| 61 | + |
| 62 | + # list_fix: z_exp is a design parameter for KvExpPiecewise, always excluded |
| 63 | + fix_set = Set(list_fix) |
| 64 | + ps.hydraulic.kv isa KvExpPiecewise && push!(fix_set, :z_exp) |
| 65 | + for (i, name) in enumerate(params.name) |
| 66 | + mask[i] && name ∈ fix_set && (mask[i] = false) |
| 67 | + end |
| 68 | + |
| 69 | + # list_sameLayer: keep only first occurrence per name |
| 70 | + same_set = Set(list_sameLayer) |
| 71 | + seen = Set{Symbol}() |
| 72 | + for (i, name) in enumerate(params.name) |
| 73 | + mask[i] || continue |
| 74 | + name in same_set || continue |
| 75 | + name in seen ? (mask[i] = false) : push!(seen, name) |
| 76 | + end |
| 77 | + params[mask, :] |
| 78 | +end |
| 79 | + |
| 80 | + |
| 81 | +# list_sameLayer: broadcast the source-layer value to all hydraulic profile layers |
| 82 | +# after update!, and before update_hydraulic! (VG: m depends on n) |
| 83 | +# list_fix: excluded from params by get_params, so update! never touches them |
| 84 | +function update_params!(ps::SoilColumn{FT,N}, paths, theta; |
| 85 | + params=nothing, |
| 86 | + list_sameLayer::Vector{Symbol}=Symbol[], |
| 87 | + list_fix::Vector{Symbol}=Symbol[]) where {FT<:Real,N} |
| 88 | + |
| 89 | + isnothing(params) && (params = optim_params(ps, :hydraulic; list_sameLayer, list_fix)) |
| 90 | + length(paths) == length(theta) || |
| 91 | + throw(ArgumentError("paths and theta must have the same length, got $(length(paths)) and $(length(theta)).")) |
| 92 | + |
| 93 | + update!(ps, paths, theta; params) |
| 94 | + |
| 95 | + # broadcast list_sameLayer from source layer to all layers before update_hydraulic! |
| 96 | + # 要求在 |
| 97 | + for (i, name) in enumerate(params.name) |
| 98 | + name in list_sameLayer || continue |
| 99 | + I = params.path[i][end] # layer index |
| 100 | + I isa Integer || continue |
| 101 | + h = getproperty(ps.hydraulic.profile, name) |
| 102 | + fill!(h, h[I]) |
| 103 | + end |
| 104 | + |
| 105 | + # update derived fields (e.g. VanGenuchten m = 1 − 1/n); must precede _sync_ksat! |
| 106 | + update_hydraulic!(ps.hydraulic.profile) |
| 107 | + _sync_ksat!(ps.hydraulic.kv, ps.hydraulic.profile, ps.hydraulic.dz_cm) |
| 108 | + |
| 109 | + # rebuild AoS caches from updated SoA profiles |
| 110 | + @inbounds for i in 1:N |
| 111 | + ps.hydraulic.layers[i] = ps.hydraulic.profile[i] |
| 112 | + ps.thermal.layers[i] = ps.thermal.profile[i] |
| 113 | + end |
| 114 | + return nothing |
| 115 | +end |
0 commit comments