From 183593d4e60bc6469b752b2ef8105e3be6511c71 Mon Sep 17 00:00:00 2001 From: Julia Sloan Date: Tue, 10 Dec 2024 17:19:27 -0800 Subject: [PATCH] add water_conservation_test artifact --- README.md | 1 + water_conservation_test/Manifest.toml | 305 ++++++++++++++++++ water_conservation_test/Project.toml | 2 + water_conservation_test/README.md | 12 + water_conservation_test/create_artifact.jl | 20 ++ .../run_water_conservation_exp.sh | 26 ++ .../water_conservation_exp.jl | 178 ++++++++++ 7 files changed, 544 insertions(+) create mode 100644 water_conservation_test/Manifest.toml create mode 100644 water_conservation_test/Project.toml create mode 100644 water_conservation_test/README.md create mode 100644 water_conservation_test/create_artifact.jl create mode 100644 water_conservation_test/run_water_conservation_exp.sh create mode 100644 water_conservation_test/water_conservation_exp.jl diff --git a/README.md b/README.md index 59ad374..108bb4d 100644 --- a/README.md +++ b/README.md @@ -60,6 +60,7 @@ S. Gupta et al 2022, 2024](https://github.com/CliMA/ClimaArtifacts/tree/main/soi - [ERA5 Monthly Averages on Single Levels 1979-2024](https:////github.com/CliMA/ClimaArtifacts/tree/main/era5_monthly_averages_single_level_1979_2024) - [Lehmann et al. 2008 Fig. 8 evaporation data](https:////github.com/CliMA/ClimaArtifacts/tree/main/lehmann2008_evaporation) - [Forty years of ERA5 forcing data for ClimaLand](https:////github.com/CliMA/ClimaArtifacts/tree/main/forty_yrs_era5_land_forcing_data) +- [Water conservation test data for RichardsModel](https:////github.com/CliMA/ClimaArtifacts/tree/main/water_conservation_test) ### Coupler/shared diff --git a/water_conservation_test/Manifest.toml b/water_conservation_test/Manifest.toml new file mode 100644 index 0000000..5ac0b9c --- /dev/null +++ b/water_conservation_test/Manifest.toml @@ -0,0 +1,305 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.11.2" +manifest_format = "2.0" +project_hash = "74753091693fecdd4a33be59e19da0ff479810d2" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.2" + +[[deps.ArtifactUtils]] +deps = ["Downloads", "Git", "HTTP", "Pkg", "ProgressLogging", "SHA", "TOML", "gh_cli_jll"] +git-tree-sha1 = "f7c3ca6c70ea6b035effe5428eb05231ba890c2d" +uuid = "8b73e784-e7d8-4ea5-973d-377fed4e3bce" +version = "0.2.4" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +version = "1.11.0" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +version = "1.11.0" + +[[deps.BitFlags]] +git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d" +uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35" +version = "0.1.9" + +[[deps.ClimaArtifactsHelper]] +deps = ["ArtifactUtils", "Downloads", "Pkg", "REPL", "SHA"] +path = "../ClimaArtifactsHelper.jl" +uuid = "6ffa2572-8378-4377-82eb-ea11db28b255" +version = "0.0.1" + +[[deps.CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "bce6804e5e6044c6daab27bb533d1295e4a2e759" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.6" + +[[deps.ConcurrentUtilities]] +deps = ["Serialization", "Sockets"] +git-tree-sha1 = "ea32b83ca4fefa1768dc84e504cc0a94fb1ab8d1" +uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb" +version = "2.4.2" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +version = "1.11.0" + +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.ExceptionUnwrapping]] +deps = ["Test"] +git-tree-sha1 = "d36f682e590a83d63d1c7dbd287573764682d12a" +uuid = "460bff9d-24e4-43bc-9d9f-a8973cb893f4" +version = "0.1.11" + +[[deps.Expat_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "e51db81749b0777b2147fbe7b783ee79045b8e99" +uuid = "2e619515-83b5-522b-bb60-26c02a35a201" +version = "2.6.4+1" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" +version = "1.11.0" + +[[deps.Git]] +deps = ["Git_jll"] +git-tree-sha1 = "04eff47b1354d702c3a85e8ab23d539bb7d5957e" +uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2" +version = "1.3.1" + +[[deps.Git_jll]] +deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"] +git-tree-sha1 = "399f4a308c804b446ae4c91eeafadb2fe2c54ff9" +uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb" +version = "2.47.1+0" + +[[deps.HTTP]] +deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] +git-tree-sha1 = "6c22309e9a356ac1ebc5c8a217045f9bae6f8d9a" +uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" +version = "1.10.13" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +version = "1.11.0" + +[[deps.JLLWrappers]] +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "be3dc50a92e5a386872a493a10050136d4703f9b" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.6.1" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.4" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "8.6.0+0" + +[[deps.LibGit2]] +deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" +version = "1.11.0" + +[[deps.LibGit2_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] +uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" +version = "1.7.2+0" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.11.0+1" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +version = "1.11.0" + +[[deps.Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "61dfdba58e585066d8bce214c5a51eaa0539f269" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.17.0+1" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" +version = "1.11.0" + +[[deps.LoggingExtras]] +deps = ["Dates", "Logging"] +git-tree-sha1 = "f02b56007b064fbfddb4c9cd60161b6dd0f40df3" +uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" +version = "1.1.0" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" +version = "1.11.0" + +[[deps.MbedTLS]] +deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"] +git-tree-sha1 = "c067a280ddc25f196b5e7df3877c6b226d390aaf" +uuid = "739be429-bea8-5141-9913-cc70e7f3736d" +version = "1.1.9" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.6+0" + +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2023.12.12" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" + +[[deps.OpenSSL]] +deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"] +git-tree-sha1 = "38cb508d080d21dc1128f7fb04f20387ed4c0af4" +uuid = "4d8831e6-92b7-49fb-bdf8-b643e874388c" +version = "1.4.3" + +[[deps.OpenSSL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "7493f61f55a6cce7325f197443aa80d32554ba10" +uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +version = "3.0.15+1" + +[[deps.PCRE2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" +version = "10.42.0+1" + +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.11.0" +weakdeps = ["REPL"] + + [deps.Pkg.extensions] + REPLExt = "REPL" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.2.1" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.4.3" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +version = "1.11.0" + +[[deps.ProgressLogging]] +deps = ["Logging", "SHA", "UUIDs"] +git-tree-sha1 = "80d919dee55b9c50e8d9e2da5eeafff3fe58b539" +uuid = "33c8b6b6-d38a-422a-b730-caa89a2f386c" +version = "0.1.4" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "StyledStrings", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" +version = "1.11.0" + +[[deps.Random]] +deps = ["SHA"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +version = "1.11.0" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +version = "1.11.0" + +[[deps.SimpleBufferStream]] +git-tree-sha1 = "f305871d2f381d21527c770d4788c06c097c9bc1" +uuid = "777ac1f9-54b0-4bf8-805c-2214025038e7" +version = "1.2.0" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" +version = "1.11.0" + +[[deps.StyledStrings]] +uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b" +version = "1.11.0" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.0" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +version = "1.11.0" + +[[deps.TranscodingStreams]] +git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.11.3" + +[[deps.URIs]] +git-tree-sha1 = "67db6cc7b3821e19ebe75791a9dd19c9b1188f2b" +uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" +version = "1.5.1" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +version = "1.11.0" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" +version = "1.11.0" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.13+1" + +[[deps.gh_cli_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "b54ff233a6a183ad366d29fa0e3ff9bb921692db" +uuid = "5d31d589-30fb-542f-b82d-10325e863e38" +version = "2.35.0+0" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.59.0+0" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+2" diff --git a/water_conservation_test/Project.toml b/water_conservation_test/Project.toml new file mode 100644 index 0000000..8e4179f --- /dev/null +++ b/water_conservation_test/Project.toml @@ -0,0 +1,2 @@ +[deps] +ClimaArtifactsHelper = "6ffa2572-8378-4377-82eb-ea11db28b255" diff --git a/water_conservation_test/README.md b/water_conservation_test/README.md new file mode 100644 index 0000000..91f7c83 --- /dev/null +++ b/water_conservation_test/README.md @@ -0,0 +1,12 @@ +# Water conservation test dataset + +This dataset is used as a reference solution for a particular setup of the +ClimaLand.jl soil water model (`RichardsModel`). It was created by running +the `water_conservation.jl` experiment with an explicit solver and a very +small timestep. To recreate it, we need to check out a commit of ClimaLand.jl +where this type of model was stepped fully explicitly; it has since been changed +to be stepped with a mixed implicit/explicit algorithm. + +Note that if the physics implemented in `RichardsModel` changes, this dataset +will need to be regenerated so it can still be used as a comparison for the +current code. diff --git a/water_conservation_test/create_artifact.jl b/water_conservation_test/create_artifact.jl new file mode 100644 index 0000000..f30981b --- /dev/null +++ b/water_conservation_test/create_artifact.jl @@ -0,0 +1,20 @@ +using ClimaArtifactsHelper + +# Set the output directory +output_dir = basename(@__DIR__) * "_artifact" +if isdir(output_dir) + @warn "$output_dir already exists. Content will end up in the artifact and may be overwritten." + @warn "Abort this calculation, unless you know what you are doing." +else + mkdir(output_dir) +end + +# Run the bash script to create the artifact +run(`bash run_water_conservation_exp.sh $output_dir`) + +@info "Data file generated!" + +# create_artifact_guided( +# output_dir; +# artifact_name = basename(@__DIR__), +# ) diff --git a/water_conservation_test/run_water_conservation_exp.sh b/water_conservation_test/run_water_conservation_exp.sh new file mode 100644 index 0000000..1cb3ddd --- /dev/null +++ b/water_conservation_test/run_water_conservation_exp.sh @@ -0,0 +1,26 @@ +#!/bin/bash + +# This script checks out a specific commit of the ClimaLand.jl GitHub +# repository, and uses it to run the provided `water_conservation_exp.jl` +# experiment. That file will run the experiments with both flux and Dirichlet +# boundary conditions, and save the solution results. +# The output data files are saved in the paths water_conserv_exp/ref_soln_flux.csv +# and water_conserv_exp/ref_soln_dirichlet.csv, relative to the directory +# where this script is run. + +# Inform the user that the script requires one input argument, if it is not provided +[[ $# -lt 1 ]] && echo "Error! Usage: $0 " && exit 0 + +# Parse the input argument for the output directory +output_dir=$1 +mkdir -p $output_dir + +# Clone the ClimaLand.jl repository +git clone https://github.com/CliMA/ClimaLand.jl.git +cd ClimaLand.jl + +# Checkout the commit that we know has what we need +git -c advice.detachedHEAD=false checkout d7cc27a9b33227cfc1aa3745a94f3ffcc595fa5c + +julia --project=experiments -e 'using Pkg; Pkg.instantiate()' +julia --project=experiments ../water_conservation_exp.jl $output_dir diff --git a/water_conservation_test/water_conservation_exp.jl b/water_conservation_test/water_conservation_exp.jl new file mode 100644 index 0000000..fe3b618 --- /dev/null +++ b/water_conservation_test/water_conservation_exp.jl @@ -0,0 +1,178 @@ +using ClimaCore +import ClimaComms +ClimaComms.@import_required_backends +import SciMLBase +import ClimaTimeSteppers as CTS +using Plots +using Statistics +using DelimitedFiles + +using ClimaLand +using ClimaLand.Soil +using ClimaLand.Domains: Column +import ClimaUtilities.OutputPathGenerator: generate_output_path + +rmse(v1, v2) = sqrt(mean((v1 .- v2) .^ 2)) + +# Define simulation times +t_start = Float64(0) +dt = Float64(0.001) +t_end = Float64(1e6) + +FT = Float64 + +stepper = ClimaLSM.RK4() +ode_algo = CTS.ExplicitAlgorithm(stepper) + +# van Genuchten parameters for clay (from Bonan 2019 supplemental program 8.2) +ν = FT(0.495) +K_sat = FT(0.0443 / 3600 / 100) # m/s +vg_n = FT(1.43) +vg_α = FT(0.026 * 100) # inverse meters +θ_r = FT(0.124) +S_s = FT(1e-3) #inverse meters +hcm = vanGenuchten{FT}(; α = vg_α, n = vg_n) + +zmax = FT(0) +zmin = FT(-0.5) +nelems = 50 + +params = Soil.RichardsParameters(; + ν = ν, + hydrology_cm = hcm, + K_sat = K_sat, + S_s = S_s, + θ_r = θ_r, +) +soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems) +sources = () + +# Set flux boundary conditions (used for calculating mass balance) +flux_in = FT(-1e-7) +top_bc = Soil.WaterFluxBC((p, t) -> flux_in) +flux_out = FT(0) +bot_bc = Soil.WaterFluxBC((p, t) -> flux_out) + +boundary_fluxes = (; top = top_bc, bottom = bot_bc) + +soil = Soil.RichardsModel{FT}(; + parameters = params, + domain = soil_domain, + boundary_conditions = boundary_fluxes, + sources = sources, +) + +exp_tendency! = make_exp_tendency(soil) +set_initial_cache! = make_set_initial_cache(soil) +imp_tendency! = make_imp_tendency(soil) +jacobian! = make_jacobian(soil) + +Y, p, coords = initialize(soil) +@. Y.soil.ϑ_l = FT(0.24) +set_initial_cache!(p, Y, FT(0.0)) + +jac_kwargs = + (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = jacobian!) + +prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction(T_exp! = exp_tendency!), + Y, + (t_start, t_end), + p, +) +sol = SciMLBase.solve(prob, ode_algo; dt = dt, saveat = dt) + +# Save flux BC reference solution artifact +savedir = ARGS[1] +soln_file = joinpath(savedir, "ref_soln_flux.csv") +open((soln_file), "w") do io + writedlm(io, parent(sol.u[end].soil.ϑ_l), ',') +end + +# TODO remove this after running once + # Calculate water mass balance over entire simulation + mass_end = sum(sol.u[end].soil.ϑ_l) + mass_start = sum(sol.u[1].soil.ϑ_l) + t_sim = sol.t[end] - sol.t[1] + # Flux changes water content every timestep (assumes constant flux_in, flux_out) + mass_change_exp = -(flux_in - flux_out) * t_sim + mass_change_actual = mass_end - mass_start + relerr = abs(mass_change_actual - mass_change_exp) / mass_change_exp + @show relerr + @assert relerr < FT(1e-9) + mass_errors[i] = relerr + + +# Perform simulation with Dirichlet boundary conditions +top_state_bc = Soil.MoistureStateBC((p, t) -> ν - 1e-3) +flux_out = FT(0) +bot_flux_bc = Soil.WaterFluxBC((p, t) -> flux_out) +boundary_conds = (; top = top_state_bc, bottom = bot_flux_bc) + +soil_dirichlet = Soil.RichardsModel{FT}(; + parameters = params, + domain = soil_domain, + boundary_conditions = boundary_conds, + sources = sources, +) + +exp_tendency! = make_exp_tendency(soil_dirichlet) +set_initial_cache! = make_set_initial_cache(soil_dirichlet) +imp_tendency! = make_imp_tendency(soil_dirichlet) +jacobian! = make_jacobian(soil_dirichlet) +update_aux! = make_update_aux(soil_dirichlet) + +rmses_dirichlet = Array{FT}(undef, length(dts)) +mass_errors_dirichlet = Array{FT}(undef, length(dts)) + +Y, p, coords = initialize(soil_dirichlet) +@. Y.soil.ϑ_l = FT(0.24) +set_initial_cache!(p, Y, FT(0.0)) + +jac_kwargs = + (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = jacobian!) + +prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction(T_exp! = exp_tendency!), + Y, + (t_start, t_end), + p, +) +# TODO remove saveat, sv after running once +saveat = Array(t_start:dt:t_end) +sv = (; + t = Array{Int64}(undef, length(saveat)), + saveval = Array{NamedTuple}(undef, length(saveat)), +) +cb = ClimaLand.NonInterpSavingCallback(sv, saveat) +sol = SciMLBase.solve( + prob, + ode_algo; + dt = dt, + callback = cb, + adaptive = false, + saveat = saveat, +) + +# TODO remove this after running once + # Calculate water mass balance over entire simulation + # Because we use Backward Euler, compute fluxes at times[2:end] + flux_in_sim = + [parent(sv.saveval[k].soil.top_bc)[1] for k in 2:length(sv.saveval)] + + + mass_end = sum(sol.u[end].soil.ϑ_l) + mass_start = sum(sol.u[1].soil.ϑ_l) + t_sim = sol.t[end] - sol.t[1] + mass_change_exp = -(sum(flux_in_sim) * dt - flux_out * t_sim) + mass_change_actual = mass_end - mass_start + relerr = abs(mass_change_actual - mass_change_exp) / mass_change_exp + @assert relerr < 1e11 * eps(FT) + mass_errors_dirichlet[i] = relerr + +# Save Dirichlet BC reference solution artifact +soln_file_dirichlet = joinpath(savedir, "ref_soln_dirichlet.csv") +open((soln_file_dirichlet), "w") do io + writedlm(io, parent(sol.u[end].soil.ϑ_l), ',') + writedlm(io, parent(sol.u[end].soil.ϑ_l), ',') +end