Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Streamline isotope code: a) use C instead of δ as state varible, b) reduce allocations, c) possibly include δ update in RHS f #80

Open
wants to merge 2 commits into
base: develop
Choose a base branch
from
Open

Streamline isotope code: a) use C instead of δ as state varible, b) reduce allocations, c) possibly include δ update in RHS f #80

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dc0ab4a
refactor: Comment out more unused code in LWFBrook90R_updateIsotopes_…
Feb 8, 2024
7c6b5e1
test: fix regression test on Mac
Feb 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/src/index.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ For further details read through the section [User Guide](@ref) or refer to sect


## Citing LWFBrook90.jl
When using LWFBrook90.jl please cite <!-- [Bernhard et al. (2020)](http://dx.doi.org/10.5194/egusphere-egu2020-17975) -->
When using LWFBrook90.jl please cite <!-- [Bernhard et al. (2024)](http://dx.doi.org/10.5194/egusphere-egu2020-17975) -->
>TODO: update citation with published article.


Expand Down
60 changes: 42 additions & 18 deletions docs/src/user-guide.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ The structure of the input CSV's is illustrated by the example input data sets `

For convenience, input CSV files can be generated from a script that sets up a simulation with the R package [LWFBrook90R (v0.4.3)](https://github.com/pschmidtwalter/LWFBrook90R#usage). Instead of running the simulation with `run_LWFB90()`, the same arguments can be used to generate the input files for LWFBrook90.jl using the R function provided in the file `generate_LWFBrook90jl_Input.R`. Note that the input file `meteoiso.csv` needs to be generated separately and the files containing the initial conditions (`initial_conditions.csv` and `soil_discretization.csv`) also need to be extended manually with the isotope values (see structure of these input files below).

To load input data and prepare a simulation follow the instructions in section [Example Script 01](@ref) or alternatively use the sample script `main_with_isotopes.jl`. NOTE: these will be replaced with Jupyter-notebooks generated with Literate.jl (TODO).
To load input data and prepare a simulation follow the instructions in section [Example Script 01](@ref) or alternatively use the sample script `main_with_isotopes.jl`. NOTE: these will be replaced with Jupyter-notebooks generated with Literate.jl.

In case you're unfamiliar with Julia, there are various ways to run a script such as `main.jl`: One possibility is to open the Julia REPL and run the script using `include(“main.jl”)`. Alternatively, the editor VS Code in combination with the Julia extension ([julia-vscode.org](https://www.julia-vscode.org)), provides a complete IDE for programming in Julia.

Expand Down Expand Up @@ -153,23 +153,27 @@ Roadmap to include calibration data intends include of:

as well as

- isotopic composition of soil water (δ_soil)
- isotopic composition of xylem water (δ_xylem)
- Isotopic composition of soil water (δ_soil)
- Isotopic composition of xylem water (δ_xylem)

Below the example structure of data sets for
`psi.csv` contains the soil matric potential (ψ):

| dates | depth_m | psi_kPa |
| ---------- | ------- | ------- |
| YYYY-MM-DD | m | kPa |
| 2019-12-23 | 0.00 | -0.1 |
| 2019-12-23 | 0.20 | -0.93 |
| 2019-12-23 | 0.40 | -0.1 |
| 2019-12-23 | 0.80 | -0.27 |
| 2019-12-23 | 1.60 | -0.1 |
| 2019-12-24 | 0.00 | -0.1 |
| ... | ... | ... |
| | | |
`psi.csv` contains the soil matric potential (ψ) of different sensor series or as site average:

| dates | depth_cm | depth_nominal_cm | psi_kPa | series |
| ---------- | -------- | ---------------- | ------- | ------- |
| YYYY-MM-DD | cm | cm | kPa | SITEAVG |
| 2022-08-30 | 15 | 15 | -5.77 | SITEAVG |
| 2022-08-30 | 200 | 150 | -2.87 | SITEAVG |
| 2022-08-30 | 300 | 150 | -0.10 | SITEAVG |
| 2022-08-30 | 50 | 50 | -1.98 | SITEAVG |
| 2022-08-31 | 80 | 80 | -0.13 | SITEAVG |
| 2022-08-31 | 15 | 15 | -7.53 | SITEAVG |
| 2022-08-31 | 200 | 150 | -2.95 | SITEAVG |
| 2022-08-31 | 300 | 150 | -0.10 | SITEAVG |
| 2021-06-17 | 50 | 50 | -3.00 | SITEAVG |
| 2021-06-17 | 80 | 80 | -0.17 | SITEAVG |
| ... | ... | ... | ... | ... |
| | | | | |

`theta.csv` contains the soil moisture (volumetric water content, θ):

Expand All @@ -186,8 +190,28 @@ Below the example structure of data sets for

`delta_soil.csv` contains the isotopic signature of soil water:

TODO...
| dates | depth_m | delta18O_permil | delta2h_permil |
| ---------- | ------- | ----------- | ------- |
| YYYY-MM-DD | cm | permil | permil |
| 2020-07-07 | 0 | -5.89 | -35.96 |
| 2020-07-07 | 15 | -8.75 | -57.99 |
| 2020-07-07 | 50 | -8.81 | -59.12 |
| 2020-07-07 | 80 | -9.72 | -65.41 |
| 2020-07-21 | 0 | -5.83 | -29.54 |
| 2020-07-21 | 15 | -8.73 | -53.58 |
| 2020-07-21 | 50 | -9.31 | -58.47 |
| 2020-07-21 | 80 | -10.14 | -64.71 |
| ... | ... | ... | ... |

`delta_xylem.csv` contains the isotopic signature of xylem water:

TODO...
| dates | species | treeID | delta18O_permil | delta2H_permil |
| ---------- | ------- | ------ | --------------- | -------------- |
| YYYY-MM-DD | - | - | permil |permil |
| 2021-08-03 | Beech | 271 | -8.416 | -72.562 |
| 2021-08-03 | Beech | 3518 | -8.946 | -76.543 |
| 2021-08-03 | Beech | 3519 | -8.484 | -72.523 |
| 2021-08-17 | Beech | 271 | -8.378 | -76.044 |
| 2021-08-17 | Beech | 3518 | -8.769 | -76.036 |
| 2021-08-17 | Beech | 3519 | -8.135 | -70.475 |
| ... | ... | ... | ... |... |
58 changes: 28 additions & 30 deletions src/func_DiffEq_definition_cb.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -796,19 +796,19 @@ function LWFBrook90R_updateIsotopes_GWAT_SWAT_AdvecDiff!(u, t, integrator)
D_¹⁸O_ᵏ⁺¹ .= D⁰_¹⁸O .* τw .* θᵏ⁺¹ .+ p_DISPERSIVITY .* abs.(q) # m²/day
D_²H_ᵏ⁺¹ .= D⁰_²H .* τw .* θᵏ⁺¹ .+ p_DISPERSIVITY .* abs.(q) # m²/day

# Define concentrations of source/sink terms in transport equation (TRANI, DSFLI, INFLI, SLVP)
### Define δ signature of in- and outflows (TRANI, DSFLI, INFLI)
# TODO(bernhard) for debugging: use p_δ18O_PREC instead of p_fu_δ18O_SLFL
δ18O_INFLI = p_δ18O_PREC(integrator.t)#TODO(bernhard): debug remove workaround and set again = p_fu_δ18O_SLFL
δ2H_INFLI = p_δ2H_PREC(integrator.t) #TODO(bernhard): debug remove workaround and set again = p_fu_δ2H_SLFL

C_¹⁸O_INFLI = LWFBrook90.ISO.δ_to_C.(δ18O_INFLI, LWFBrook90.ISO.R_VSMOW¹⁸O, LWFBrook90.ISO.Mi_¹⁸O)
# C_¹⁸O_TRANI = C_¹⁸Oᵏ # no fractionation occurring, i.e. outflux composition equal to storage composition
# C_¹⁸O_DSFL = C_¹⁸Oᵏ # no fractionation occurring, i.e. outflux composition equal to storage composition
C_²H_INFLI = LWFBrook90.ISO.δ_to_C.(δ2H_INFLI, LWFBrook90.ISO.R_VSMOW²H, LWFBrook90.ISO.Mi_²H)
# C_²H_TRANI = C_²Hᵏ # no fractionation occurring, i.e. outflux composition equal to storage composition
# C_²H_DSFL = C_²Hᵏ # no fractionation occurring, i.e. outflux composition equal to storage composition
# TODO(bernhard): for TRANI and DSFL we are instead using Cᵏ⁺¹ (i.e. use the concentrations in the LHS in the implicit scheme...)
# # Define concentrations of source/sink terms in transport equation (TRANI, DSFLI, INFLI, SLVP)
# ### Define δ signature of in- and outflows (TRANI, DSFLI, INFLI)
# # TODO(bernhard) for debugging: use p_δ18O_PREC instead of p_fu_δ18O_SLFL
# δ18O_INFLI = p_δ18O_PREC(integrator.t)#TODO(bernhard): debug remove workaround and set again = p_fu_δ18O_SLFL
# δ2H_INFLI = p_δ2H_PREC(integrator.t) #TODO(bernhard): debug remove workaround and set again = p_fu_δ2H_SLFL

# C_¹⁸O_INFLI = LWFBrook90.ISO.δ_to_C.(δ18O_INFLI, LWFBrook90.ISO.R_VSMOW¹⁸O, LWFBrook90.ISO.Mi_¹⁸O)
# # C_¹⁸O_TRANI = C_¹⁸Oᵏ # no fractionation occurring, i.e. outflux composition equal to storage composition
# # C_¹⁸O_DSFL = C_¹⁸Oᵏ # no fractionation occurring, i.e. outflux composition equal to storage composition
# C_²H_INFLI = LWFBrook90.ISO.δ_to_C.(δ2H_INFLI, LWFBrook90.ISO.R_VSMOW²H, LWFBrook90.ISO.Mi_²H)
# # C_²H_TRANI = C_²Hᵏ # no fractionation occurring, i.e. outflux composition equal to storage composition
# # C_²H_DSFL = C_²Hᵏ # no fractionation occurring, i.e. outflux composition equal to storage composition
# # TODO(bernhard): for TRANI and DSFL we are instead using Cᵏ⁺¹ (i.e. use the concentrations in the LHS in the implicit scheme...)

### Compute δ signature of evaporating flux (SLVP)
δ¹⁸O_SLVP = u_δ18O_SWATI[1] # TODO(bernhard): disabling evaporation fractionation (TODO: should yield solution similar to Stumpp 2012 model)
Expand All @@ -823,12 +823,12 @@ function LWFBrook90R_updateIsotopes_GWAT_SWAT_AdvecDiff!(u, t, integrator)
# Cᵢ_SLVP = ( (Cᵢ - ε¹⁸O_eq)/α¹⁸O_eq - h*δ¹⁸O_a - ε¹⁸O_dif ) /
# (1 - h + ε¹⁸O_dif/1000) # [‰]

C_¹⁸O_SLVP .= 0
C_²H_SLVP .= 0
C_¹⁸O_SLVP[1] = LWFBrook90.ISO.δ_to_C(δ¹⁸O_SLVP, LWFBrook90.ISO.R_VSMOW¹⁸O, LWFBrook90.ISO.Mi_¹⁸O)
C_²H_SLVP[1] = LWFBrook90.ISO.δ_to_C(δ²H_SLVP, LWFBrook90.ISO.R_VSMOW²H, LWFBrook90.ISO.Mi_²H)
# E¹⁸O = C_¹⁸O_SLVP * aux_du_SLVP[1] * 0.001 # kg/m3 * mm/day * 0.001 m/mm # (kg/m²/day¹)
# E²H = C_²H_SLVP * aux_du_SLVP[1] * 0.001 # kg/m3 * mm/day * 0.001 m/mm # (kg/m²/day¹)
# C_¹⁸O_SLVP .= 0
# C_²H_SLVP .= 0
# C_¹⁸O_SLVP[1] = LWFBrook90.ISO.δ_to_C(δ¹⁸O_SLVP, LWFBrook90.ISO.R_VSMOW¹⁸O, LWFBrook90.ISO.Mi_¹⁸O)
# C_²H_SLVP[1] = LWFBrook90.ISO.δ_to_C(δ²H_SLVP, LWFBrook90.ISO.R_VSMOW²H, LWFBrook90.ISO.Mi_²H)
# # E¹⁸O = C_¹⁸O_SLVP * aux_du_SLVP[1] * 0.001 # kg/m3 * mm/day * 0.001 m/mm # (kg/m²/day¹)
# # E²H = C_²H_SLVP * aux_du_SLVP[1] * 0.001 # kg/m3 * mm/day * 0.001 m/mm # (kg/m²/day¹)

### Prepare terms to evaluate linear system to be solved
Δt = integrator.t - integrator.tprev # days
Expand Down Expand Up @@ -856,13 +856,13 @@ function LWFBrook90R_updateIsotopes_GWAT_SWAT_AdvecDiff!(u, t, integrator)

# # Setup linear system (implicit solver) (From Braud et al. 2005)

# Define boundary conditions
### TOP: surface isotopic flux E = 0 (as evaporation in LWFBrook is implemented as sink term (see above E¹⁸O))
BCFlux_top¹⁸O = 0 # kg/m²/day¹
BCFlux_top²H = 0 # kg/m²/day¹
### BOTTOM:
BCFlux_bottom¹⁸O = q[N] * C_¹⁸Oᵏ[N] # kg/m²/day¹
BCFlux_bottom²H = q[N] * C_²Hᵏ[N] # kg/m²/day¹
# # Define boundary conditions
# ### TOP: surface isotopic flux E = 0 (as evaporation in LWFBrook is implemented as sink term (see above E¹⁸O))
# BCFlux_top¹⁸O = 0 # kg/m²/day¹
# BCFlux_top²H = 0 # kg/m²/day¹
# ### BOTTOM:
# BCFlux_bottom¹⁸O = q[N] * C_¹⁸Oᵏ[N] # kg/m²/day¹
# BCFlux_bottom²H = q[N] * C_²Hᵏ[N] # kg/m²/day¹

# # Make matrix for implicit time stepping (solving Ax = b) (see Braud et al. 2005)
# # TODO(bernharf): preallocate this matrix A (and vectors...)
Expand Down Expand Up @@ -1132,8 +1132,6 @@ function LWFBrook90R_updateIsotopes_GWAT_SWAT_AdvecDiff!(u, t, integrator)
aux_du_TRANI.*Cᵢ²H_TRANI
)

#TODO(bernhard): Diffusion slows the code down considerably and makes it unstable (for NLAYER=81 only, NLAYER=41 seems okay...)

# # NOTE: below max(0.001,u_SWATIᵏ⁺¹) makes the code more robust
# du_Cᵢ¹⁸_SWATI = -C_¹⁸Oᵏ./max.(0.001,u_SWATIᵏ⁺¹) .* dVdt .+ 1 ./ max.(0.001,u_SWATIᵏ⁺¹) .* (
# -diff¹⁸O_upp*1000 .+ diff¹⁸O_low*1000 .+ qCᵢ¹⁸O_upp .- qCᵢ¹⁸O_low .+
Expand Down Expand Up @@ -1187,8 +1185,8 @@ function LWFBrook90R_updateIsotopes_GWAT_SWAT_AdvecDiff!(u, t, integrator)
du_δ2H_GWAT = LWFBrook90.ISO.dxdt_to_dδdt.(du_Cᵢ²H_GWAT, Cᵢ²H_GWAT, LWFBrook90.ISO.R_VSMOW²H)
end
# go back from atom fraction to delta values
du_δ18O_SWATI .= LWFBrook90.ISO.dxdt_to_dδdt.(du_Cᵢ¹⁸_SWATI, C_¹⁸Oᵏ, LWFBrook90.ISO.R_VSMOW¹⁸O)
du_δ2H_SWATI .= LWFBrook90.ISO.dxdt_to_dδdt.(du_Cᵢ²H_SWATI, C_²Hᵏ, LWFBrook90.ISO.R_VSMOW²H)
du_δ18O_SWATI .= LWFBrook90.ISO.dxdt_to_dδdt(du_Cᵢ¹⁸_SWATI, C_¹⁸Oᵏ, LWFBrook90.ISO.R_VSMOW¹⁸O)
du_δ2H_SWATI .= LWFBrook90.ISO.dxdt_to_dδdt(du_Cᵢ²H_SWATI, C_²Hᵏ, LWFBrook90.ISO.R_VSMOW²H)

# 5) Apply changes explicitly
# Since we update this in the callback we can't overwrite `du` and let DiffEq.jl do
Expand Down
2 changes: 1 addition & 1 deletion test/03-regression-tests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -480,7 +480,7 @@ end
fname_illustrations = "out/$git_status_string/TESTSET_DAV2020modified-regressionTest-FLUXES"
mkpath(dirname(fname_illustrations))

pl_fluxes = plot_simulated_fluxes_vs_reference(simulated_fluxes, loaded, d_out);
pl_fluxes = plot_simulated_fluxes_vs_reference(df_simulatedFluxes, loadeddf);
Plots.plot!(legend = :topright, size=(2000,1000), layout = (7,5))
savefig(Plots.plot(pl_fluxes, size=(2000,1400), layout = (7,5), dpi=300), fname_illustrations*"_fluxes.png")
end
Expand Down
Loading