Skip to content

Commit

Permalink
feat: API replace get_δ() with get_states()/get_fluxes()
Browse files Browse the repository at this point in the history
  • Loading branch information
Fabian Bernhard committed Apr 5, 2024
1 parent bc68d5d commit ef9cecf
Show file tree
Hide file tree
Showing 4 changed files with 10 additions and 25 deletions.
5 changes: 3 additions & 2 deletions examples/example-script-01.jl
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,8 @@ using LWFBrook90
## - `simulate!()`
## - `get_soil_()`
## - `get_states()`
## - `get_δ()`,
## - `get_fluxes()`,
## - `get_forcing()`,

# Define simulation model by reading in system definition and input data from input files.
# When printed, the generated SPAC model gives a summary.
Expand Down Expand Up @@ -210,7 +211,7 @@ show(df_out_daily)
## if true # simulate_isotopes
## df_δsoil = get_soil_([:δ18O, :δ2H],
## example_result; depths_to_read_out_mm = depth_to_read_out_mm)
## δ_results = get_δ(simulation_mod)
## δ_results = get_states(simulation_mod)
## end
##
## timepoints = simulation_mod.ODESolution_datetime
Expand Down
1 change: 0 additions & 1 deletion src/LWFBrook90.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@ export RelativeDaysFloat2DateTime
# read out results for soil domain variables
export get_soil_
# read out results for aboveground/scalar variables
export get_δ, get_delta # get_mm #
export get_forcing, get_states, get_fluxes # TODO(bernhard): make sure we have documentation for these exported variables
# read out results for water partitioning
export get_water_partitioning
Expand Down
22 changes: 2 additions & 20 deletions src/func_postprocess.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ function get_states(simulation::DiscretizedSPAC; days_to_read_out_d = nothing) #
timepoints =
isnothing(days_to_read_out_d) ? unique(round.(simulation.ODESolution.t)) : # case nothing: return daily by default
days_to_read_out_d == :integrator_step ? simulation.ODESolution.t : # if requested return for each integration step stored in ODESolution
(1:10)[1] isa Number ? days_to_read_out_d : # else assume we got a vector of days
days_to_read_out_d[1] isa Number ? days_to_read_out_d : # else assume we got a vector of days
error("Unknown `days_to_read_out_d` provided.")

# i) scalar states (+PREC forcing)
Expand Down Expand Up @@ -132,7 +132,7 @@ function get_fluxes(simulation::DiscretizedSPAC; days_to_read_out_d = nothing) #
isnothing(days_to_read_out_d) ? unique(round.(simulation.ODESolution.t)) : # case nothing: return daily by default
days_to_read_out_d == :integrator_step ? error("Cumulative fluxes should not be read out on subdaily intervals.") : # if requested return for each integration step stored in ODESolution
# days_to_read_out_d == :integrator_step ? simulation.ODESolution.t : # if requested return for each integration step stored in ODESolution
# (1:10)[1] isa Number ? days_to_read_out_d : # else assume we got a vector of days
days_to_read_out_d[1] isa Number ? days_to_read_out_d : # else assume we got a vector of days
error("Unknown `days_to_read_out_d` provided.")
t_ref = simulation.ODESolution.prob.p.REFERENCE_DATE

Expand Down Expand Up @@ -443,24 +443,6 @@ end

##########################
# Functions to get values either isotopes or amounts (or otherwise)
"""
get_δ(simulation::DiscretizedSPAC; days_to_read_out_d = nothing)
Returns a DataFrame with the isotopoic compositions of the inputs and state variables: PREC, GWAT, INTS, INTR, SNOW, RWU, XYLEM.
By default, the values are returned for each simulation timestep.
The user can define timesteps as `days_to_read_out_d` by optionally providing a numeric vector, e.g. saveat = 1:1.0:100
"""
function get_δ(simulation::DiscretizedSPAC; days_to_read_out_d = nothing)
compartments_to_extract = [:PREC, :GWAT, :INTS, :INTR, :SNOW, :RWU, :XYLEM]
units_to_extract = [:d18O, :d18O, :d18O, :d18O, :d18O, :d18O, :d18O]
dfd18O = get_scalars(compartments_to_extract, units_to_extract, simulation, days_to_read_out_d);

units_to_extract = [:d2H, :d2H, :d2H, :d2H, :d2H, :d2H, :d2H]
dfd2H = get_scalars(compartments_to_extract, units_to_extract, simulation, days_to_read_out_d);

return hcat(dfd18O, dfd2H[:, Not(:time)])
end
get_delta = get_δ

function get_scalars(compartments_to_extract, units_to_extract, simulation::DiscretizedSPAC, days_to_read_out_d = nothing)
solution = simulation.ODESolution;
Expand Down
7 changes: 5 additions & 2 deletions test/03-regression-tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,10 @@ function get_some_states_to_compare(example_result)
# reduce(hcat, [example_result.ODESolution[t_idx].accum for t_idx = eachindex(example_result.ODESolution)])
)
u_mm = hcat(DataFrame(time = t_out), DataFrame(permutedims(mat_aboveground), df_aboveground_names[:]))
u_δ = get_δ(example_result; days_to_read_out_d = t_out)
u_δ = innerjoin(
get_states(example_result, days_to_read_out_d = t_out),
get_fluxes(example_result, days_to_read_out_d = t_out),
on = [:dates, :PREC_d2H, :PREC_d18O])[:,[:dates,:PREC_d18O,:GWAT_d18O,:INTS_d18O,:INTR_d18O,:SNOW_d18O,:RWU_d18O,:XYLEM_d18O,:PREC_d2H,:GWAT_d2H,:INTS_d2H,:INTR_d2H,:SNOW_d2H,:RWU_d2H,:XYLEM_d2H]]
## vector quantities
u_belowground = get_soil_([:SWATI, :W, , , :K, :δ18O, :δ2H], example_result;
days_to_read_out_d = t_out)
Expand Down Expand Up @@ -167,7 +170,7 @@ function test_states_comparison(u_mm, u_δ, u_belowground, currSPAC,
# Test scalar states
compare_scalar = (A,B; nans = false) -> all(isapprox.(Matrix(A), Matrix(B); nans))
@test compare_scalar(loaded_u_ref, u_ref)
@test compare_scalar(loaded_u_δ[:, Not(:time)], u_δ[:, Not(:time)]; nans=true)
@test compare_scalar(loaded_u_δ[:, Not(1)], u_δ[:, Not(1)]; nans=true) # Not(1) removes :time or :dates

# Test vector states
compare_vector = (A,B) -> all(isapprox.(Matrix(A), Matrix(B)))
Expand Down

0 comments on commit ef9cecf

Please sign in to comment.